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10 . ABSTRACi 


The National Aeronautics and Space Administration is sponsoring a two-phase program to 
evaluate the capability of current computational fluid dynamics codes for use as a de- 
sign and analysis tool in gas turbine engine combustor development, and to improve the 
quantitative accuracy of such codes as necessary. 

This report describes the Phase I study to evaluate the physical modeling embodied in 
the codes. The objectives were to identify shortcomings in the models and to provide a 
program plan to improve the quantitative accuracy. The physical models studied were 
for: turbulent mass and momentum transport, heat release, liquid fuel spray, and gas- 
eous radiation. The approach adopted was to test the models against appropriate bench- 
mark-quality test cases from experiments in the literature for the constituent flows 
that together make up the combustor real flow. 

The results of the Phase I effort indicate that considerable work remains to be done 
before the computational fluid dynamics codes can be recognized as viable design and 
analysis tools. This work covers both physical models and numerical methods. 
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FOREWORD 


This document was prepared by the Conmercial Engineering Division of the Pratt 
& Whitney Aircraft Group of United Technologies Corporation, East Hartford, 
Connecticut, to describe an evaluation of the suitability of contemporary 
computational fluid dynamics codes to serve as design and analysis tools in 
the development of con^ustors for gas turbine engines. The evaluation was 
accomplished under the National Aeronautics and Space Administration-sponsored 
Aerothermal Modeling Program, Phase I, Contract NAS3-23524. The work was 
performed under the direction of National Aeronautics and Space Administration 
project Monitor, Dr. Edward J. Mularz. 

Key Commercial Engineering Division contributors to this program were: Dr. 
Geoffrey J. Sturgess, Technical Program Manager and Principal Investigator, 
and Dr. Saadat A. Syed and Mr. Keith R. McManus, both of Powerplant Analysis. 

Subcontract support in assembling the experimental data base was furnished by 
United Technologies Research Center, a Division of United Technologies 
Corporation; key United Technologies Research Center contact was Dr. John B. 
McVey. Numerous colleagues of Dr. McVey contributed their services as 
gatekeepers in the areas of their technical expertise. Mr. Louis M. Chiappetta 
also computeriiod the data base. 

This report has been assigned the Contractor's report number PWA-5907-19. 
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1,0 SUMMARY 


The National Aeronautics and Space Administration Is sponsoring a two-phase 
program to evaluate the capability of current computational fluid dynamics 
codes for use as a design and analysis tool In gas turbine engine combustor 
development, and to Improve the quantitative accuracy of such codes as 
necessary. 

This report describes the Phase I study to evaluate the physical modeling 
embodied In the codes. The objectives were to Identify shortcomings In the 
models and to provide a program plan to 1n^)rove the quantitative accuracy. 

The Pratt & Whitney Aircraft 2D-TEACH and 3D-TEACH computer codes have been 
used as the basis of the study, representing as they do, state-of-the-art 
modeling. 

In Phase 1 the physical models studied were for: turbulent mass and momentum 
transport, heat release, liquid fuel spray, and gaseous radiation. The 
approach adopted was to test the models against appropriate benchmark-quality 
test cases from experiments In the literature for the constituent flows that 
together make up the combustor real flow. 

The results of the Phase I effort Indicate that considerable work remains to 
be done before the computational fluid dynamics codes can be recognized as 
viable design and analysis tools. This work covers both physical models and 
numerical methods. 


2.0 INTRODUCTION 


Conputatlonal fluid dynamics codes have shown potential for application to the 
Internal flows represented by the gas turbine engine combustion chamber 
(References 1 and 2). The motivations for such an application are powerful 
(Reference 3). One man's reasonably balanced view of the place of such 
numerical modeling In the engineering of a combustor has been presented by 
Odgers (Reference 4)» and an approach to modeling combustors has been 
described (Reference 5). However, although calculations of specific combustion 
chamber flows using particular models have been made, e.g.. References 6 and 
7, and many comparisons with measurements exist for flows relevant to 
combustors, e.g.. References 8 and 9, no complete, systematic study has been 
earned out concerning the performance of these models. 

In response to the recommendations of References 1 ..nd 2, the Acrothermal 
Modeling Program is being undertaken. The overall objective of the program is 
to provide as a design tool a quantitative capability in aerothermal modeling 
for gas turbine engine combustion chambers. 

The objectives of this Phase I study are to provide a complete assessment of 
the available numerical modeling capability, and to recommend a comprehensive 
program for model improvement. The model studied and the assessment of It are 
reported here. The reewnmended program for model improvement is reported 
separately. 
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3.0 PROGRAM DESCRIPTION 


As part of a large research project entitled Turbine Engine Hot Section 
Technology (HOST), the Aerothermal Modeling Program was initiated by the 
National Aeronautics and Space Administration at the Lewis Research Center to 
provide a quantitatively accurate numerical modeling capability for the design 
and development of gas turbine engine combustion chambers. This effort is 
planned in two phases, consisting of an assessment of modeling capability and 
improvement of the accuracy of modeling. The Phase I study, which is described 
in this report, consisted of two major tasks: 

Task 1 - Model and Benchmark Definition, involved describing the models to 
be evaluated, establishing a data base of relevant experiments, 
and selecting benchmark quality test cases from the data base. 

Task 2 - Model Assessment, where the models were run for the test cases 
and the calculations compared with the measurements, assessment 
of the resulting performance including identification of 
deficiencies, and, a program plan was prepared for model 
improvement . 


3.1 Program Goals 

Phase I specific program goals are as follows: 

1. To define a detailed computational fluid dynamics code and describe 
its constituent modules and their interaction to show the usefulness 
of the output. 

2. To define an assessment procedure and assemble a structured 
experimental or analytical data base that is compatible with this 
procedure, identifying shortcomings in the data base. 

3. To select from the data base experiments of benchmark quality against 
which the model can be compared. 

4. To assess the model calculations against the experimental 
measurements to identify weaknesses. 

5. To formulate a comprehensive program plan that seeks to improve the 
model to quantitative accuracy. 


3 


3.2 Program Schedule 

The Phase I program covered a period of 9 months. Of this effort, 4 months 
elapsed time was df'^oted to assembling the data base, and 5 1/2 months to 
assessing the modei >erformance. There was swne overlap of these tasks. The 
period of performance and scope of the study were such that an extremely 
detailed analysis of the specifics of any particular model failure could not 
be carried out. 


3.3 Program Approach 

The program was directed towards state-of-the-art numerical modeling and was 
to critique all relevant physical modeling, numerical methods and grid 
systems. The experimental data base was to be formed from either experiments 
in the literature or, experiments performed on a real gas turbine combustor. 
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4.0 MODEL DESCRIPTION 


A brief description of the models is given. This description covers the 
modeling of the equations, the physical models that provide closure to the 
equations, the discretization of the equations in finite difference form, and 
the solution algorithms. Most of the developments are presented in 
two- dimens ions in the interest of brevity; extension to three-dimensions 
follows logically. 

The model selected for this study is derived from the TEACH codes developed 
initially as a teaching aid at Imperial College, London, England, (Reference 
10). The acronym TEACH (Teaching Elliptic Axisymmetric Characteristics 
Heuristically) represents a generic solution technique.H’his model was used in 
Bie study as two computer programs developed by Pratt & Whitney Aircraft: 
2D-TEACH and 3D-TEACH. The 2D-TEACH code can solve axisymmetric or planar 
geometries, while 3D-TEACH can solve fully three-dimensional problems. These 
codes represent current production state-of-the-art calculations. Relation to 
the design system is given. 

4.1 Outline of Calculation 

Figure 1 presents a flow diagram of the calculation process. It is divided 
into four major sections: assembly and modeling of equations, physical 

modeling, computer solution, and output. 

The currently most practical approach to solving a complex turbulent flow is 
to stay within the framework of continuum mechanics and to use a statistical 
description of the turbulence, coupled with the accepted Eulerian flow 
description provided by the Navier-Stokes equations of motion. Hence, an 
instantaneous quantity in the Navier-Stokes equations is described as the sum 
of a time-averaged value and a randomly fluctuating value. 

When the statistical description of an instantaneous quantity is substituted 
into the Navier-Stokes equations and time-averaged, the resulting equation set 
is known as the Reynolds equations (Reference 11). These equations are similar 
to the Navier-Stokes equations except that time-averaged quantities are used, 
and for the appearance of time-averaged correlations of fluctuating quantities. 

Turbulent motions increase the apparent viscosity of a fluid by some orders of 
magnitude. If laminar diffusion terms are therefore very much smaller than 
turbulent diffusion terms, then neglect of fluctuations in laminar viscosity 
is acceptable, it is a frequently used practice (Reference 12) to also neglect 
terms involving fluctuating density, although this implies that temperature 
differences in the flow are not large. For variable density flows, an 
alternative, density-weighted decomposition of the conservation equations 
(Favre averaging) is sometimes used. The choice of weighted or unweighted 
equations is presently somewhat arbitrary (Reference 13). Favre-averaging 
requires modeling terms like p'u’ which then appear in the equations. Since 
this presents some difficulties, unweighted equations were used. 
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Figure 1 Flow Diagram of Calculation Process 
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The simplifications of neglecting fluctuating laminar viscosity and 
fluctuating density terms give Reynolds equations that are expressed in terms 
of ti me-mea n quantities and cross-correlations of fluctuating velocities such 
as P u' -ju' j. These are the Reynolds stresses, and they result in a closure 
problem. Tdrbulence modeling provides the necessary descriptions of the 

Reynolds stresses in known or determinable quantities. 

When the flow consists of more than one chemica l species, modeling is also 
required for the turbulent mass fluxes?^ u'jm'|. These terms arise from 
applying the statistical treatment of turbulence to an instantaneous species 
transport equation. The instantaneous energy equation is given the same 
treatment. 

If the flow is two-phase, the presence of the liquid or solid phase can be 
accounted for by introducing a void fraction term, ft. This allows coupling of 
the phases since when ft is unity the equations revert to their original form 
and when 0<ft<l the gaseous phase responds to the presence of the other phase. 

The modeled equations can be algebraically manipulated into a general form 
that simplifies solution. As an example, for steady state flow in 
two-dimensional cylindrical coordinates: 


0x 




+ — s- (rp ftv $) 
rdr ' 


A. 

8x 




+ 


3 

r3r 






( 1 ) 


where, 

^eff,<^ 




^d,d> 


any of the independent variables 

an appropriate turbulent exchange coefficient, depending on 
what represents 

a so-called "source term" for the gaseous phase which lumps 
together all other terms in a given equation not included in 
the first four terms in Equation 1. 

a droplet source term for a liquid phase, and, 


subscript 'g' denotes gas phase 
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Figure 2 gives the appropriate set of source terms for each <f>. To obtain these 
terms the chemical species transport equation was applied to a simple, 
irreversible, one-step chemical reaction. Closure was provided by a 
two-equation turbulence model, requiring transport equations for the kinetic 
energy of turbulence K, and its dissipation E, and the assumption that 
turbulent heat and mass diffusivities can be related to mOTjentum diffusivity 
by turbulent Prandtl and Schmidt numbers, a*. In Figure 2, Sfy denotes the 
time-averaged reaction rate of the fuel, and Qj^ is the net volumetric heat 
exchange due to radiation. 

To complete the equation set an equation of state is provided: 


P = 


R T 


( 


"•ox 

Mw 


ox 


"’fu 

m 


m, 


PR 


fu 


Mw 


PR 


) 


where Mw is the molecular weight. 


( 2 ) 


The temperature term is calculated in an iterative manner from the composition 
and enthalpy values using temperature versus enthalpy tables for the 
individual species. 


4.2 Physical Models 

Modeling is used to describe the relevant physical processes in order to 
provide closure to the equations. A turbulence model is required to provide 
the turbulent momentum and mass fluxes, the chemical reaction rate is found 
from a combustion model, a spray model yields droplet source terms, and the 
heat transfer due to radiation is given by a radiation model. 

The models used in the present versions of the TEACH series of codes are 
defined and described briefly below. The particular models presented do not 
represent the only choices. They were selected on a pragmatic basis with the 
demands of the intended application to the fore, and with deference to the 
state of the modeling. 


4.2.1 Turbulence 


In order to solve the Reynolds equations for the mean quantities the Reynolds 
stresses have to be described. The so-called two equation, or K-E, turbulence 
model is currently used. 

Reynolds stresses can consist of two parts - a shear stress and a normal 
stress - on the faces of an elemental volume in the flow. 
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To address the shear stress portion of the Reynolds stress Boussinesq's 
analogy is used. In laminar flow. Stoke' s law relates the shear stress In a 
flowing fluid to the viscosity of the fluid through the velocity gradient. 
Turbulence increases the apparent viscosity of the fluid by some orders of 
magnitude due to transfer of energy from the mean flow through the eddy 
cascade to the smallest eddies, where it is dissipated eventually on the 
molecular level as heat. Boussinesq drew an analogy with laminar flow to 
relate shear stress in a turbulent fluid to an "eddy viscosity", Mt* With 
eddy viscosity modeling, an expression is then needed for 


The Prandtl -Kolmogorov definition of eddy viscosity (References 14 and 15) is 
used currently. Since specific (eddy viscosity per unit mass) has to have 
dimensions of L^/T, it can be written. 


l*^ (3) 

where, Vj = velocity scale related to the turbulence 
1 = length scale appropriate to Vj, and, 

m, n = unknown indices. 

If the specific kinetic energy of turbulence is defined as 


K = 2^u'^ + v'^ + w'^) (4) 

then, Vj equal VT is appropriate, and 1 then characterizes the size of the 
eddies containing the turbulence energy, where, 

1 = 1^ ; Ij = C/4 K^/^/e 


and, C /4 = constant 

e = dissipation rate of turbulence 
Hence, with the regular definition of 

/Xj = P ^ (5) 

The normal strains on an elemental volume due to turbulence pressure can be 
represented by summing the net fluctuating specific dynamic pressures in each 
direction, i.e.. 


P 


t 


1 

I 





) 
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Using Equation 4, this reduces to, 
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far 

Q^'ality 


= 2/3 K (6) 

Hence, the Reynolds stresses can be written In two-dimensions, using tensor 
notation. 


- p u' .u' . = 
^ 1 j 



2/3 p K 8. . 

ij 


(7) 


where 8^j Is the Kronecker delta such that 1 equals j denotes a normal 
stress and 8 equals unity, while 1 not equal to j denotes tangential stress 

only and 8 equals zero. 

For variable density flows an additional term of -2/3Pt6-j jCut (ap/3xi)]/P 
appears on the right hand side of Equation 7. 

Equation 5 supplies 

In Equation 5 the eddy viscosity has been expressed in terms of K and s. 
Transport equations for these can be derived as follows: For K, multiply each 
of the instantaneous momentum equations by its respective fluctuating velocity 
component, sum the three resulting equations, and then time-average. For e, 
differentiate the instantaneous momentum equations, multiply throughout by 
P3ui/3xi, sum, and time-average, (Reference 16). Some additional modeling 
of fluctuating terms is then required. 

The modeled form of the K-equatlon 1s: 



Similarly the form of the e-equatlon Is: 


Convection Diffusion 


Production 


3pu .e 



+ 






u. 

1 


u . 


J 



(9) 


Dissipation 

2 

" r c 
' ^ K 
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The model contains five empirical constants, C^, Cc], and Cc2- 

"Standard" values for these constants are given in Table 1. Zfj. comes from 

near-wall shear stress and turbulence energy measurements, where production 
and dissipation of turbulence are nearly equal; it is consistent with a law of 
the wall (log-law). Cc2 comes from the decay of grid turbulence, Cg] comes 
from near-wall turbulence considerations, and and cr^ come from 
"computer optimization" to give good answers over a range of flows. 


Table 1 

Standard Values for Turbulence Model Constants 


Constant 

Value 

Cfi 

0.09 

Ccl 

1.44 

Cc2 

1.92 


1.0 


1.3 


For variable density flows the turbulent heat and mass fluxes must be 
calculated. The flowfield is found from an effective turbulent eddy viscosity. 
It is consistent to also base the turbulent heat and mass transfer on an 
effective thermal or mass eddy diffusivity. 

The eddy diffusivity approach gives for the flux of a scalar. 


( 10 ) 


where 0is a scalar such as temperature or specie concentration. This is the 
so-called gradient hypothesis where the turbulent transport of ^ is 
proportional to the mean gradient of 6, 

The turbulent eddy diffusivity (or the exchange coefficient F) is physically 
always the same order of magnitude as the turbulent kinematic eddy viscosity 
and their ratio is called the turbulent Prandtl or Schmidt number 
where , 




( 11 ) 
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In wean field closures the isotropic turbulent Prandtl/Schwidt number has to 
be specified. Knowledge of used as a means of finding ^rom an 
already calculated value of obtained using Equation 5. Thus, transport of 
momentum is taken as the basic gradient hypothesis transport process, rather 
than transport of a scalar quantity. 

A recommended value for has been 0.9. This is now considered appropriate 
for boundary layer flows only, and a value of 0.7 has been recommended for 
general flows (Reference 13). A value of 0.5 has been recommended for 
recirculation zones, (Reference 17). 


4.2,2 Fuel Spray 

The fuel spray model has to provide information on the droplet source terms in 
Equation 1. This information is required for each cell formed by the finite 
difference mesh upon which the equation is solved. 

The approach used is to calculate droplet trajectories through the flow field. 
This method involves embedding the droplet equations of motion in Lagrangian 
form into the Eulerian framework of the gas-phase equations. 

The spray is represented by individual droplets, termed the "computational 
droplet." Each computational droplet represents a "parcel" of like droplets 
all having the same initial size, velocity vector and temperature. All the 
droplets in the actual spray are represented by a size distribution. A 
Rosin-Ranmler distribution is computationally convenient to use. The actual 
spray size distribution is sampled statistically to produce the computational 
droplets; real fuel flow rate is maintained at its correct value. The 
computational droplets are fired into the gas field. 

The droplet equation of motion is due to Soo (Reference 18) with his 
subsequent simplifying reductions, together with the assumption that the fuel 
vapor produced by the droplet has an initial velocity equal to that of the 
droplet. This yields. 


dUd . 

"d- Tf- - "d ^<“9.1 - “d,i> " V1 


where, 


p-j = Q _!_9_ 
d d 


(u_ - u. .) 

g,i d,i 


( 12 ) 


(13) 


i 


Ud i = droplet velocity in x, r and 6 directions, respectively, i.e., 

Ud, Vd, Wd 

rd = droplet radius (droplets assumed spherical) 

Cq = drag coefficient (relationship of Williams, Reference 19) 
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Then, the droplet position at any Instant of time is given by. 


dx 

dt 




( 14 ) 


The turbulent diffusion of droplets is accounted for by a stochastic approach 
(Reference 20). In the droplet equation of motion. Equation 12, is 
interpreted as the instantaneous value of gas velocity, i.e. * 

^9>i ^ ^g.i ^'g.i 

The time-mean value Uq -j is obtained from the flow field solution. The 
instantaneous value Ug can -be obtained using random sampling of an assumed 
Gaussian probability attribution function for turbulence energy, and the 
turbulence is taken as being isotropic. The standard deviation of Uq ^ is 
then, 


a = 



(16) 


K being the specific kinetic energy of turbulence, which is available 
throughout the field from the gas-phase solutions. The fluctuating velocity 
u'g,i can then be obtained from Equation 5. Hence, with isotropy, 

Ug = u + u' ; Vg = V + u* ; Wg = w + u' 

The time interval over which the droplet equations are integrated is the time 
which a particular instantaneous velocity interacts with the droplet. To 
determine this time interval it is assumed that an instantaneous velocity is 
associated with a particular turbulent eddy. The droplet interacts with the 
eddy until either the droplet traverses the eddy or, the eddy is broken up, 
whichever occurs first. The length scale of the eddy is obtained from the 
turbulence model as. 


3 

= Cfi K ^ /S 

and the transit time is determined analytically from a simplified and 
linearized form of the droplet equation of motion. 

The droplet source terms that arise from this Particle Source in Cell (PSIC) 
technique (Reference 21), involve for each computational cell. 
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(i) mass transfer between droplets and gas (continuity equation), 
( 11 ) momentum transfer (momentum equations), 

(ill) enthalpy exchange (energy equation), 

( 1 v) specie exchange (specie equations), 

(v) radiation exchange (energy equation) , 

Expressions for these source terms may be written as follows: 


S 


d,m 


I 1 



( 17 ) 


where, liidi, 111^2 
M 
N 


= mass flux of liquid droplets at inlet and outlet of cell 
= number of droplets representing size range L 
= total number of size ranges representing the spray 


and, 


m - II ^ n3 i 
"’dl ■ Z ^dl \ 

( 18 ) 

*"d2 " ^ ^d ^d2 "l 


where, \ 
Dd 

Simi larly, 



number of real droplets per second In size range L crossing cell 
diameter of droplet 


N 


E 




L =1 


I 

H 



K = 1 


['”> 0,2 

['”> 0,2 


'”>0,1] K 
<”>o,i]k 


L 


L 


( 19 ) 
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d,w 


K=1 


M j. 

B Z- I'^’d.Z ■ <“>d,lj 

K=1 “* 


( 21 ) 


’d.h 


’d,ml 


where. 

Finally, 


N 

Z 

L=1 




K=1 


^"’d,! ■ ‘ ^evap^ K 


m) 


L 

I 

l“i 


= Z 


M S ^*'l"’dl ■ '"l"’d2^ K 
K=1 


(23) 


fraction of specie in droplet, (F] « 1.0 for pure fuel and 
zero for oxidant and products, normally; see combustion model). 


N 


4a 


’d,R “ T 


L=1 


M 


E (T ^ - T 
‘'d ^'d 'g - 


B Z <-1 ‘l r/> 

K=1 ^ K 

X. 

where, t^ = time required for the droplets to cross cell 
Eji = emissivity of the droplets 

T(j = temperature of droplets 

Tg = gas temperature 

a = Stefan - Boltzman constant 

V = cell volume 


(24) 


Phenomenological modeling is used to provide closure to Equations 17 and 
19-24. Provision is made for fractional distillation of a multi-component fuel. 
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At a wall the droplet can be either reflected, or remain on the wall. If it 
remains on the wall it is taken out of the calculation. If it is reflected, 
new initial velocity conditions are specified. This means that the sign of the 
velocity component normal to the wall is altered, and tracking is resumed. A 
droplet is only ’^'•"lected when it crosses an axis of symmetry. Note that the 
analysis takes no account of either droplet collision or secondary breakup 
after the spray is formed. 

4.2.3 Combustion 

The time-averaged chemical reaction rate Sfu in Figure 2 is provided by the 
combustion model. 

The current model is very single. It neglects molecular transport and chemical 
kinetics completely, and has been classified as a hydrodynamics only 
procedure. The argument is that in turbulent flames the time scale of chemical 
reaction is much smaller than the time scale of turbulence and, therefore, 
reaction rate is determined by mixing of fuel and air since for all practical 
purposes, the chemical reaction rate is infinitely fast. Hence, rate of 
reaction is influenced only by the time scale of turbulence. This is the 
familiar "mixed is burned" hypothesis, and is appropriate for the gas turbine 
combustor when operating at high power conditions. 

The time-average reaction rate appearing in the source terms is obtained from 
the model of Magnussen and Hjc'-tger (Reference 22). This approach deals with 
the inhomogeneities which appear in the mixture as a consequence of the 
interactions between turbulence and chemical reaction, through a modification 
of Mason and Spalding's eddy breakup model (Reference 23). 

The modified eddy breakup model as applied presently can be described as 
follows. In lean diffusion flames fuel and oxidant occur in separate eddies, 
and the rate of reaction is determined by the rate of eddy dissipation. 
Consequently, there will be a relationship between turbulent fluctuations in 
the flow and the mean concentration of active species. Accordingly, the 
volumetric reaction rate of fuel can be expressed, 

= PAi mf^J (e/K) ' (25) 

where, ifu "" local time-mean mass fraction of fuel 

e/K = reciprocal of eddy life time 

Ai = constant 

i 

In regions of the flame where the time-mean mixture is rich, the oxidant will 
be the reacting species that shows the greatest intermittency. Hence, in this 
case it is the dissipation of the oxidant-bearing eddies that limit reaction 
rate. Hence, 
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m 


Yu 


PA 


ox / 

1 ITlKi 


where, = local time-mean mass fraction of oxidant. 

In premixed flames eddies will contain both fuel and oxidant simultaneously, 
and such eddies will be separated by eddies containing hot products. Reaction 
rate is determined in this case by flame spread, and by dissipation of 
hot-eddies when the concentration of hot products is low. These two processes 
were combined into a single equation of the same form as Equations 25 and 26, 


— — / \ 

Vu " ^'^ i'^2 n+^ ' (f ) 

where, mpy. = local time-mean mass fraction of products 
A2 = constant 

In general, even diffusion flames will contain some level of premixing, so 
Equations 25-27 will apply simultaneously. Thus, the reaction will proceed at 
a rate determined by the lowest of the three rates. 

The values used for the constants are those recommended by Magnussen and 
MJertager, i.e.. 


A] = 4.0 ; A]A2 = 2.0 - A3 


(28) 


4.2.4 Radiation 

The radiation model provides the heat transfer due to radiation, Q|^. The 
approach adopted is a relatively simple one, consistent with the simple 
canbustion model. 

It is assumed that the radiating medium is grey. The radiation transport 
equation for this medium then becomes independent of the frequency of 
radiation (Reference 24). This equation when integrated over the 47T solid 
angle yields the following "energy equation" for the radiation field, 

= CK (Em-Er) (29) 

where 

Er, = Radiation energy density = -2 / I dco 

^ y 47T 
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C = speed of l^ght 

O' =' Stefan-Boltzman constant 

T = temperature 

I = radiation intensity which a function of direction and position 

o) = solid angle 

K = Absorption coefficient 

n = direction vector 

If the radiation transport equation is multiplied by the direction vector and 
integrated over Att solid angle, the following "momentum equation" for the 
radiation field is obtained. 


where Pp 


CV • Pr = - K Fr 


Radiation Pressure 



Ion* do> 


(30) 


To close the above set of equations (Equations 29, 30} ^ is expressed in 
terms of by expanding the radiation intensity in a Taylor series in n 
and truncating it after the first two terms. This yields a relation between 
radiation energy density and radiation pressure which when substituted in 
Equation 30 yields 




Putting Equation 31 in Equation 29 gives 

This equation, in the presence of absorbing/emitting particles, can 
in cylindrical coordinates as: 

4 4 

8 1 ^R . 1 8 r 0^R . . 4 . 4\ p . 4 4\ . 


(31) 


be written 
0 (32) 
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4 / C 

Tj^ = 1. 0/4(7 y Idw = ^ X radiation energy density. 

I = Radiation intensity 

w = Wave length 

<7 = Stefan-Boltzman Constant 

Tg = Gas temperature 

Tp = Particle temperature 

K = Absorption Coefficient 

Ke = Effective absorption Coefficient, Ke = K + Kp + 4/3 Dh 
Dh = Hydraulic diameter of the enclosure 

Ep = Particle emissivitv 

The wall boundary condition can be written as: 


T ^ 

1 a'R _ 1 

■JiT an " ■? 
e 



(33) 


where 

Zyi - Emissivity of the wall 
= Reflectivity of the wall 
n = Direction normal to the wall surface 
= Wall temperature 

The solution of the modeled radiation transport equation. Equation 32, 
requires the specification of the gas absorption coefficient K. The evaluation 
of K for combustion gases over the range of temperatures and pressures 
encountered in the combustor can be very complex. The gas absorption 
coefficient is therefore calculated as a frequency averaged absorption 
coefficient. It can be shown, that 


dE. 


K = 


m 


dx 


X = 0 


where En, = Emissivity of the gas mixture. 


(34) 
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The problem now reduces to calculating the emissivlty of the combustion gases. 

Barteld's (Reference 25) model Is Incorporated. This model Is simple but 
general enough so that constants used to curve-fit data for one set of 
operating conditions can be used under a different set. 


This model assumes that the combustion gases are made up of a mixture of four 
gray gases whose emissivlty can be calculated as: 


4 

^m ■ E ®g,n 
n»l 




(35) 


where 


X » path length, length of the control volume 


g.n 




Kp = constant 

p = partial pressure of products of combustion 

The constants and b]], etc., are given In Table 2. These constants were 
determined by making a curve-fit for the combustion products of methane. The 
fit was made over the temperature range 1800®R (1000“K) to 3600®R (2000°K) and 
a partial pressure - path length range of 0.01 - 4.0 ft. atniospheres (0.003 - 
1.2 bar meter) - at an average pressure of one atmosphere. This model should 
be adequate for predicting radiation from gaseous fuel fired combustors at low 
pressures. 


This model was modified to take Into account collision broadening due to high 
ambient pressures. The following modification was made: 


= 
m 


4 

= I. 

n=l 


g.n 




('■ 




1/4 


px 


) 


(36) 


where P is the ambient pressure. The constants remain unchanged. 

From Equations 34 and 36 the absorption coefficient needed for the radiation 
transport equation can be obtained as: 




4 


z 

n= 1 





(37) 


21 


Table 2 


OF POUR QIIkuITY 

Constants of Four Gray Gas Model for the Products of Methane Combustion 


n 


1 

2 

3 

4 


‘>ln 

h 

‘’2n* 

^3n 

+ .218 

+ .283 

+ 

+ .535 

-.240 

+ 

+ .131 

+ .089 

- 

+ .116 

-.132 

+ 


10® 

b 

®4n- 


.065 

-.026 

0.106 

.009 

+ .002 

1.653 

.127 

+ .032 

21.212 

.053 

-.008 

146.308 


4.3 Solution Procedure 

An outline of the solution procedure for the set of equations (Figure 2) 
represented by the general transport equation of Equation 1 is provided in 
sufficient detail that appreciation can be gained of how the elements of this 
procedure might influence solution accuracy. 

Rearrangement of the equations into the general form represented by Equation 1 
enables one solution algorithm to be used for all equations. The equation set 
is solved using a steady state, implicit, finite difference numerical 
procedure. An initial guess is made of the field variables, and these guesses 
are iteratively updated until the solutions have converged. Convergence is 
deemed to have been obtained when the absolute sum of the residuals over the 
whole grid of each variable goes below a specified value. 

The solution procedure is the standard TEACH approach: A hybrid 
(upwind/central) finite differencing scheme (Reference 26) is used to 
discretize the equations. Velocities and pressures are obtained through the 
SIMPLE algorithm (Reference 27) and the finite difference equations are solved 
by the ADI method (Reference 28) using the TDMA. 

4.3.1 Discretization of the Equations 

The finite difference analog of the differential equations is obtained by 
overlaying a computational mesh on the flow domain to be calculated, and 
obtaining the basic finite difference form of the partial derivatives for 
every node of the mesh from a control volume approach, (Reference 29). The 
finite difference expressions, when substituted back into the differential 
equations, yield a set of linearized, algebraic equations for every node of 
the mesh. 

There are a number of ways to derive the finite difference analogs of partial 
differential equations. The control volume or finite volume approach, because 
it is based on the satisfaction of macroscopic physical laws such as 
conservation of mass, momentum and energy, usually leads to more accurate 
results. Figure 3 illustrates the mesh and the control volume established 
about a considered mode, P. The conservation property is essential when 
combustion is taking place. 
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Figure 3 2D- Control Volume for the Finite Difference Scheme 


The generalized transport equation for variable 0 has a source cerm $ 0 , (see 
Equation 1 and Figure 2, for example). This term is linearized thus, 


^0 = A 0p + B 

and integrated over the control volume as, 

J{h 0p + B) dV = A 0p • vol + B • vol = Sp 0p + 
where. 


(38) 




vol = ydV = J y r dr dx for two- dimens ions in cylindrical coordinates. 
V xw rs 


The remainder of Equation 1 is also integrated over the control volume and 
added to Equation 38 to give, 

Cj ■ V - “w % - V * “n <^n ■ V 


(39) 


- D- (0 - 0 ) + (S 0 + S ) 

s ^s' ' p^p u' 
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where, C^, C^, etc., are convection coefficients, ^ 

Cf = ( Pu)eae ; = ( Pvlpap etc. 

ag, ap being cell face areas, and 

D£, D^, etc., are diffusion coefficients. 



When certain weighting factors are introduced. Equation 39 can be manipulated 
with the aid of continuity, and normalized, to give the form, 

Ap </>|^ + + A^^ (40) 

where, 

Ap = A|^ + A$ + Ay + Ae - Sp (41) 


and, A|y|, A 5 , etc., are functions of the convection and diffusion 
coefficients and the weighting factors. 

Equation 40 is the finite difference form of Equation 1 for <^. To solve it 
evaluation of the weighting factors is required. The difference scheme used to 
evaluate the weighting factors is based on the hybrid differencing scheme of 
Spalding (Reference 26). 


The hybrid differencing scheme is unconditionally stable and the solution is 
bounded. It uses second order central differencing^for convection and 
diffusion fluxes when the absolute value of cell Peclet number is less than or 
equal to two. When Peclet number is greater than two, first order upwind 
differencing is used for convection fluxes, and diffusion fluxes are neglected 
altogether. The switch of differencing is done both locally and directionally 
in the computational grid. Peclet number defines the relative importance of 
convective and diffusive transport, and. 


Peg = C£/D£ ; Pey, = Cy/Dy, etc. 
Applying the scheme it can be shown that for example, for 

Pg<2 ; A^ = Cn 


"2 < Pg < 2 ; A|\| 


Dm - Cy/2 


Pg >2 ; A|\| = 0 


(42) 
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Similarly* for the other A's. 

Pe<2 ; As 
-2 < Pe<2 ; As 
P g > 2 ; As 

Pe<2 ; Ae 
-2 < Pg<2 ; Ae 
P g>2 ; Ae 

Pg < 2 ; Ay 




Ds + Cs/2 

cs 


De - Ce/2 
0 


0 


"2 < Pg<2 ; Ay = Dy + 
Pg>2 ; Ay = Cy 


(43) 


(44) 


(45) 


The spatial differencing of the convective terms of the conservation equations 
in an Eularian coordinate system can result in numerical diffusion occurring. 
Use of a higher order differencing scheme eliminates or significantly reduces 
this diffusion. However, use of central differencing introduces an oscillatory 
behavior into the solution. This "wiggling" can lead to nonphysical behavior 
(Reference 30). The use of an upwind or donor-cell technique eliminates 
wiggling; however, this is accomplished by the introduction of a 
diffusive-1 ike term into the difference equations. Thus, while "numerical 
damping" suppresses oscillation, it leads to significant additional diffusion 
of the convected parameter. For flows with combustion, these parameters might 
be species concentration, temperature, etc. Unfortunately, diffusion of these 
quantities is responsible in a physical sense for flame propagation. 

Therefore, a severe restriction can be placed on the quality of quantitative 
prediction (Reference 31). 


To successfully use the hybrid differencing scheme for complicated flows, care 
must be taken in establishing the computational grid upon which the 
calculations are performed. The approximations of the algebraic expressions 
used to represent the partial differential equations become asymptotically 
exact as the distance between the nodes set up by the grid, and used to link 
the algebraic expressions, is reduced. In the limit, the number of nodes can 
be increased until an asymptote to the solution to the differential equations 
is achieved. In practice, this increase is limited by computer storage and the 
cost of the calculation. However, it is not just the number of nodes that are 
used which is important in determining the accuracy of a solution, but also 
the distribution of those nodes within the flow field to be determined 
(References 13, 32). This nodal distribution is important because whenever 
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curvature of the flow in the streamwise coordinate direction exists, a 
truncation error arises in the solution (Reference 33). In addition, there is 
also a problem in multidimensional flows of streamline-to-grid skewness 
(Reference 34). With upwind differencing, these effects start to have a 
damaging effect on solution accuracy when the P6cl6t number exceeds two. 

4.3.2 Pressure Algorithm SIMPLE 

The acronym SIMPLE stands for Semi-I^mplicit Method for Pressure Linked 
Equations. This algorithm is formulated to obtain a pressure field and correct 
the velocity field in such a manner that continuity and momentum equations are 
simultaneously satisfied. The need for a pressure algorithm arises because 
pressure does not appear as an independent variable in any of the transport 
equations. 

The procedure adopted is the following: The pressure field, P*, is guessed 

and the momentum equations are solved using this field. The resulting velocity 
field u* contains the right vorticity but may or may not satisfy continuity. 
Hence, an equation for the potential function, P', is derived whose solution 
yields velocity and pressure corrections, u^' and P', respectively, which 
preserve vorticity and satisfy continuity. The derivation of this equation 

requires a relationship between u^' and P' to be found. The following 
relationship is assumed. 


pu. 


= - r. 


ap‘ 

3x. 


(46) 


where is deduced from the linearized momentum equations. 
The continuity equation can then be written as 


3(pu^.) . 9p(u.* + u. ' ) 


9x, 


8x. 


9pu.* 


ax. 


3P' 

ax. 


Hence 


where D* = 


3pu*. 




i 


3 r 

3x. M 3x^. 


- D* 


(47) 


(48) 

(49) 


Equation 47, which is known as the pressure correction equation, is cast in 
the form of a generalized transport equation (the convection term is set to 
zero), and is solved usirrg the algorithm which is used to solve all the other 
transport equations. Knowing P', u^' is obtained from Equation 46. 
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The set of finite difference equations, Equation 40, are solved by means of 
the following line by line iteration method. This method known as the 
Alternating Direction I_mplicit Method (Reference 28) was initially formulated 
Tor unsteady equations; its adaptation to steady state equations is sometimes 
also known as the Alternating Direction Uerative Method. 

If the <f> values on two of the four neighboring nodes of a considered node P 
are assumed as known, then Equation 39 can be recast as 

Ap0p = A^^n + 


where 


C'*A£4>^ + Ay^y + Su (51) 

In Equation 50, neighboring values and are assumed as known. This 
assumption reduces the finite difference equation to a three point equation 
which can now be solved using any fast matrix inversion technique. Once the 
equation set is solved the direction is reversed and now and d>$ are 
assumed as known, which produces the following equation: 

Ap«p = + As<>s + C” (52) 

where 

C ' = S 5u (53) 

Equation 52 can also be solved by matrix inversion. In the above equations the 
most recent values of <f> are used. 

4.3.4 TDMA 

The ADI method reduces the finite difference equation. Equation 40, which is a 
five point equation to two, three point equations, Equations 50 and 52, which 
are solved successively using the Tri Diagonal Matrix Algorithm (TDMA). This 
algorithm is explained below. 

Equation 50 can be rearranged for the point as 

= Bj O + Cj ^ Dj 


where 

Bj = A|^/Ap > Cj = Aj/Ap 
Dj = + A^ + Sjj)/Ap 
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The points on the computation grid range from 1 to Nj in the N-S direction 
with points 1 and N: on the boundaries. Since the boundary values and 
^NJ “ known, equations for <i >2 to(^NJ-l are solved. The set of equations 
then becomes: 


^2 “ ^3 **2 + Dz 

^3 * ®3 ^4 ^ ^3 ^2 ^ ^3 ( 54 ) 


■ * • * 

•^Nj-1 “ %•-! ^Nj %-l '^Nj-2 °Nj-l 

Now since <i>] is known </>2 can be eliminated from Equation 54 and so on, 
yielding a general recurrence relation 

^j+1 °o’ 

To get the relation for Aj and Dj' Equation 55 is written as 

= Aj-i 0j + Dj-.i 


Now putting in the value of <l>j from Equation 50 


Ap As . Aj_i 


0j + ] + 


Ap - As Aj_i J j 


(55) 


(56) 


Comparing Equation 54 and 56 yields coefficients for the recurrence formula 
where Aj = A|\j/ (Ap - As Aj.]) j (57) 

Dj = (As Dj-i + Cj ) / (Ap - As Aj-i ) ^ (58) 


Using Equations 57 and 58, <^j can be calculated from Equation 55. Having 
solved for (f>^ on one N-S line, (^,-'s on the next N-S line are solved and so 
on until the'^entire solution domain is swept. The same treatment is then 
applied in the W-E direction (Equation 52). It is usually necessary to sweep 
between 1 and 3 times per iteration for optimum solution cimv->. 

4.3.5 Under-Relaxation 

Since the finite difference equations are nonlinear in nature the convergence 
is facilitated and sometimes divergence is avoided by under-relaxing the value 
of (f> being calculated as: 
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where F is an under-relaxation factor which is less than one. 

The way in whicn the above relation is introduced into the numerical procedure 
is as follows: 

sR = $u + (1 - F) Af 

It can easily be shown that the effect of introducing the above modifications 
is to under-relax according to Equation 59. 

Some skill is required to select the Mnder-relaxation factors that are best 
suited to a particular problem. 

4.3.6 Solution Procedure 

The following procedure is used to obtain the solution: 

1. Guess fields for all variables. 

2. Assemble coefficients of momentum equations and solve for U* and V* 
using prevailing pressures. 

3. Solve the pressure correction equation and update velocities and 
pressures. 

4. Solve equations for other variables. 

5. Update fluid properties such as viscosity and density. 

6. Test for convergence. If not attained use prevailing fields as new 
guesses and repeat from step 2 until convergence is attained. 

4.4 Computer Codes 

The equations to be solved, the physical models and the numerical schemes 
described above are embodied in computer programs. There are two codes 
relevant for this study. They are known as 20-TEACH and 3D-TEACH. 

The acronym TEACH is used to denote the generic solution procedure. The prefix 
2D or 3D defines the dimensional capability of the code; 2D implying planar or 
axisymmetric flows can be calculated, and 3D implying full three-dimensional 
capability. The internal organization of the two codes is similar, and the 
interactive menus, commands and cautions used are also similar so that a user 
can operate either code with equal facility. The 3D-TEACH code can regress for 
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a two-dimensional problem to produce the same solution as 2D-TEACH. Either 
Interactive or batch running Is available. These codes exist as production 
codes. This means that they are documented and maintained, and are 
user-friendly. 

Both codes are generalized such that complex geometric boundaries and Inflows 
can be represented without reprogramming. The coordinate system used Is 
orthogonal with choice of either cylindrical or cartesian coordinates. The 
basic finite difference method causes the calculation mesh to be constrained 
by the coordinate system. With a chosen coordinate system therefore, the 
geometry has to be "discretized" to fit this system. Hence, complex geometries 
have to be represented by a series of "stair-steps". The consequences of this 
are that wall shear stresses and heat transfer cannot be correctly calculated 
for curved walls or walls inclined to the major axes. Inflow and outflow 
through the elements of the stair-steps Is permitted. In addition, solid 
bodies can be placed Inside the flowfleld and may contain mass sources or 
sinks within them. Therefore, within the constraints of the coordinate system, 
the codes have great geometric flexibility. There is a choice of either 
English or S.I. units in both codes. The 3D-TEACH code is provided with cyclic 
boundary conditions for handling repeating segments (as in an annular 
combustor) . 

The planned use of the codes is for calculation of the bulk flows in the 
combustor. Adjacent to solid boundaries the local Reynolds number of the flow 
based on local velocity and distance from the wall becomes very small, and the 
turbulence model, which was developed for high Reynolds numbers, becomes 
inappropriate. For this reason and reasons of economy, the calculation does 
not proceed right through the boundary layer to the solid wall; boundary layer 
wall functions are used instead, and are patched on to the bulk solution. 

The codes are provided with appropriate pre- and post-processors for the 
convenience of the user. The pre-processor allows the user to interactively 
create or edit a flow geometry and automatically builds an input file. For 
2D-TEACH a graphics post-processor generates profile, contour, streamline and 
streakline plots; a quasi-three-dimetis’'nnal representation of profile plots to 
show flow development is also available;. This post-processor is extended for 
3D-TEACH to include a production version of Brigham Young University's 
geometric modeling system MOVIE. BYU. The streakline routine and the spray 
model are available to generate dynamic representations of the flow and 
droplet trajectories during an interactive terminal session with the codes. 

The 2D-TEACH and 3D-TEACH codes represent state-of-the-art numerical modeling 
capability, and are available as production computer programs in regular 
engineering use. 
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4.''. Relation to Design System 

A TEACH session on the computer is the computational equivalent to performing 
a physical experiment or test. All the data, and more, of an actual experiment 
or test is available for whatever subsequent manipulation is required. In 
fact, the computational experiment is infinitely more flexible than physical 
experimentation. 

In a single terminal session, the dynamics of the isothermal flow field may be 
studed using streaklines, and, using a color terminal, the experience is 
exactly like observing a water-analogy rig flow visualization. Combustion can 
then be Lurned on and the hot flow field observed in similar fashion; the 
trajectories of fuel droplets injected into and contributing to this field can 
be observed to study penetration and evaporation if desired. Profiles of 
axial, radial and tangential mean velocities, temperatures, pressures, 
densities, mass fractions of fuel, oxidant and combustion products can be 
displayed at positions throughout the field. They can also be assembled into a 
quasi-three-dimensional plot so the field development can be readily 
displayed. If desired, contour plots may be selected. Post-processors can 
generate from the field information global parameters such as outlet 
temperature pattern and profile factors, mixedness, recirculation zone size, 
strength and mass recirculated, pressure recovery, loss coefficients, etc., as 
desired. For more detailed studies turbulence quantities, and surface shear 
stress and heat transfer coefficients can be produced. 

The information generated is available immediately off the terminal screen as 
semi-permanent hardcopy, within 2-3 hours as permanent hardcopy through a 
Calcomp plotter, or as mic'^ofiche for storage. The entire output can be saved 
in a restart file for subsequent operations. 

Parametric studies of the effect of operating conditions or geometric 
variations are easily accomplished using the pre-processor to modify the 
input. There are of course, no experimental facility or instrumentation 
limitations to restrict the range of operating conditions that may be run, and 
there are no extensive delays for hardware rework. 

The literally staggering potential of computational fluid dynamics, as 
indicated above, suggests how it might be used in a design system. 

Figure 4 shows a possible revision to the design system where a module for a 
comprehensive Computational Fluid Dynamics Model is added. The orocess begins 
with the specification of design criteria and ends with certification of the 
engine. Within this design system the module labeled Computational Fluid 
Dynamics is apparent as exerting a monitoring and controlling function 
extending from preliminary layout through engine test-development and 
encompassing all stages between. Note that rig testing is now limited to 
proof-testing, and that engine development testing is much more effective due 
to the insight fed back from the added module. Note also that the system is 
now predictive: performance information becomes available at or before the 
preliminary layout stage, and diffuser and shroud flow calculations provide 
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boundary conditions to the combustor calculations. The Liner Cooling Analysis 
module is now provided with realistic flow conditions on either side of the 
combustor liner as boundary conditions allowing improved calculation of 
temperature distributions over the liner surface. 

In Figure 4 the Computational Fluid Dynamics Model module does not replace the 
Preliminary Design System module. It is presently believed that the 
combination of phenomenological modeling and mathematical modeling en^odied in 
the two modules working together represents the cost effective approach to 
design. 

Figure 5 provides an illustration of the potential use of the codes with an 
axisymmetric study of the dump region in the dump-diffuser of a combustion 
chamber. Flow is represented in this two-dimensional axisymmetric calculation 
by streamlines. The use of stair-steps is shown clearly. 
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5.0 ASSESSfCNT PHILOSOPHY 


The sub-models and the computer code embodying them have been briefly 
described above In order to define what Is to be evaluated. To carry out a 
realistic performance evaluation of these models an unambiguous approach !$ 
required. 

Real combustor flows are three-dimensional, multiphase, turbulent and 
chemically reacting flows. Detailed measurements of real ccmibustor flows are 
costly and present many severe difficulties. Analyses o^ such flows require 
that all of the effects present In the flow be treated simultaneously. All of 
the physical models - together with the numerics - are exercised. Therefore, 
it Is nearly impossible to separate effects In order to assess the performance 
of the individual elements in this situation. If the performance of individual 
models cannot be assessed, then it is impossible to Identify specific 
shortcomings and to make recommendations for In^rovements. 

It was concluded that individual models should be examined one at a time In 
progression. Experiments of Increasing complexity that require Invoking 
multiple models should be approached In sequence so that the proqresslve 
performance is evaluated. Cognizance of the effects of the numerics on 
solution accuracy would be maintained. 

It was decided that the type of experiment which Is most suitable for use In 
evaluating model elements is one conducted with high quality hardware under 
laboratory conditions where operating variables can be controlled, and where 
extensive and careful measurements are made and the experiment is documented 
fully and carefully. Complex flows should only be examined to determine the 
extent to which Individual model elements interact and to establish that the 
solution procedures are workable in complex flows. 


6.0 DATA BASE DESCRIPTION 


To evaluate the modeling a set of suitable experiments is required. Such 
experiments have been organized into a data base. 

6.1 Class of Data 

The class of experimental data to be sought and used to establish the data 
base against which the performance of the modeling can be assessed was 
determined by the decisions to eschew real combustor flows and to rely on 
laboratory quality experiments. 

The general philosophy adopted to test the models required that the flows in 
the gas turbine combustor be considered in terms of their main constituent 
processes. These processes are: 

1. Thermal radiation 

2. Fuel spray development and flow field interaction 

3. Combustion reactions and kinetics 

4. Fluid mechanics 

5. Soot formation 

Each main process was subdivided as appropriate. For example, fluid mechanics 
was broken down into the constituent flows as follows: 

a) co-axial jets, 

b) swirling jets, 

c) jets in sudden expansions, 

d) jets in crossflow, 

e) wall jets, 

f) bluff -body recirculations 

This procedure required that experiments be located in the literature for each 
of the constituent subdivisions. In addition, redundant experiments were 
required to ensure generality of any conclusions drawn about a particular 
flow/model performance. 

6.2 Assessment Procedure 

Experiments of a particular constituent flow were assessed for their 
suitability for inclusion in the data base against a number of criteria 
established to define an ideal experiment. These criteria encompassed minimum 
dimensionality, well-behaved flow, continuous variation of test parameters, 
known boundary conditions, progression in flow complexity and extensive 
instrumentation. Use of flow visualization as a diagnostic aid was also 
considered desirable. A suitable experiment might also be an analytic solution 
or an alternative but established calculation procedure. 


1) Minimum necessary flow dimensionality . Experiments in which the flows can 
be represented as onV or two-dimensional are desirable. Three-dimensional 
flow situations should be used only when specifically testing three-dimension- 
al flow modeling capability. 

2) Well-behaved flows . Flows in which instabilities, periodicity, or changes 
in gross behavior occur as flow conditions change are to be avoided. For 
example, flows in which the location of reattachment points of separated flow 
regions could undergo significant shifts as Reynolds number is changed over 
the range of interest should be avoided (unless, of course, this is the flow 
feature being tested). 

3) Continuous variation of test parameters . Experiments that are conducted 
over a wide range of values of test parameters rather than at isolated sets of 
conditions. 

4) Known boundary conditions . Entrance and exit flow profiles must be 
specified as completely as possible. Velocity, temperature and pressure 
profiles are required. Concentration profiles are important for reacting flow 
experiments; initial droplet size, velocity and spatial distributions are 
required for two-phase flows. For assessing turbulence models, initial 
profiles of turbulence intensity and integral length scale are vital. 

5) Progression in flow complexity . The ideal experiment for assessing flow 
models would consist of a series of experiments of increasing flow ccxnplexity 
such that the creditability of the analysis can be checked in stages. For 
example, the initial tests for a given flow geometry should be 
single-component, isothermal flow visualization tests which could be used to 
check the fluid mechanic aspects of the analysis. Two-component gaseous flows 
could be used to check the ability to predict mass diffusion; thermal 

diffusion could be checked using gases introduced with different initial 
conditions, reacting gaseous mixtures the next stage. A similar progression of 
experiments with two-phase flows can be constructed. 

6) Extensive instrumentation . Flow mapping experiments in which nonintrusive 
techniques are used to characterize flows throughout the chamber volume are 
highly desirable. Estimates of instrument precision should be made and 
possible sources of bias identified. Redundant measurements performed with 
different instruments are valuable. 

Another class of data is available for use in validating the computational 
accuracy of the model. These data are in reality results of computations 
carried out using well-developed and proven codes which hava been in use over 
the years to treat simple flows wherein one particular physical phenomenon is 
of interest. Results of calculations performed using tne combustor flow 
analysis given the proper boundary conditions must agree with the 
well-developed codes if the computational aspects of the model are to be 
considered valid. 


A subjective numerical rating was applied to each experiment considered 
suitable. The scale ran from 0 (poor) to 10 (good), and the score depended on 
how the experiment measured up against the criteria above. 

The survey was a limited effort. Undoubtedly there exist other works which 
might qualify as benchmark quality test cases. Experiments receiving a score 
were included in the data base, and those receiving no numerical rating were 
not included. For those experiments that received rankings, the assessments of 
how well the experiment fulfilled the various criteria are given in Table 3. 
Experiments reviewed and not included in the data base are given in Table 4. 

It should be pointed out that the individual experiments may not have been 
conducted for the purpose of generating data to be used to validate combustor 
aerothermal flow models, and therefore the degree to which the experimental 
results meet the criteria, as reflected by the numerical ranking, should not 
be construed as a general score on the quality of the experiment. 

A number of soot formation experiments were reviewed; however, it was 
concluded that these were not suitable for the present purposes and they were 
not pursued further. A new and continuing experiment by Santoro et al appears 
to hold the promise for the future (Reference 69). 
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Tabulation of Experitiients Not in Data Base 


Source 

Ref 

Category 

Subcategory 

Description 

Crabb 

41 

Isothermal 

Jets in crossflow 

Multiple Jets at 90” to 
crossflow-effects of 
non uniform Jet spacing 

Ramsey & Goldstein 

51 

Isothermal 

Jets in crossflow 

Single heated Jet at 
variable angle to 
crossflow 

Roquemore et al 

52 

Exothermic 


Propane in air diffu- 
sion flame at base of 
bluff body 

Roquemore et al 

53 

Exothermic 


Explores nonstationary 
behavior of above 

Starner & Bilger 

54 

Exothermic 


Horizontal hydrogen Jet 
flame in rectangular 
airflow with moveable 
upper and lower walls 

Kent & Bilger 

55 

Exothermic 


Continuation of above 

Starner & Bilger 

56 

Exothermic 


Extension of original 
v.ork 

Bilger S Beck 

57 

Exothermic 


Supplements Kent S 
Bilger 

Depsky 

58 

Exothermic 


Methane/air flame in 
coaxial swirling jets 

Oven, Gouldin & 
McLean 

59 

Exothermic 


Continuation of earlier 
work by Depsky 

Hutchinson et al 

44 

Exothermic 


Natural gas/air furnace 

Baker et al 

60 

Exothermic 


Continuation of 
Hutchinson et al 

Hutchinson et al 

61 

Exothermic 


Continuation of above 

Founti et al 

62 

Exothermic 

Liquid fuel 

Continuation of 
Hutchinson et al 

El Banhawy & 
Whitelaw 

63 

Exothermic 

Liquid fuel 

Rotating cup injector 
in noniwirling air- 
stream with sudden ex- 
pansion 

Brum & Samuelson 

64 

Exothermic/ 

Isothermal 


Propane/C02 into air- 
three coaxial streams 
with middle stream 
swirling 

Noyce et al 

65 

Exothermic 


Propane/air in stylized 
GT combustor 

Yule et al 

66 

Spray 


Twin fluid atomizer 
with kerosene into co- 
axial air Jets 

Modaress et al 

67 

Spray 


Glass beads into co- 
axial stream of air 

Soloman et al 

68 

Spray 


Air atomizing oil spray 
into quiescent environ- 
ment 



6.3 Organization 

The collected data were organized into the following categories: 

1) Isothermal flow - Flow field structure 

Data on confined flows with and without swirl, with and without an 
initial sudden expansion. 

2) Isothermal flow - Jet Penetration 

Data on confined flows into which penetrate single jets or rows of 
jets. 

3) Exothermic flow 

Data on confined flows with heat release. 

4) Sprays 

Data on the spatial and size distributions of sprays. 

5) Radiation 

Data on measured radiation heat fluxes in furnaces. 

An example of the collated data for a particular experiment is given in Table 
5. All of the information contained on a descriptive input sheet is contained 
in the data base, together with tables of the measured quantities. 

6.4 Management 

The collected and collated data sets of interest are contained in a data base 
stored in structured format on a magnetic tape. The data base is readily 
accessible through a user-generated retrieval system. By this means the data 
base is made readily transportable, and it can be installed on any 
moderately-sized computing system. The retrieval system is not transportable 
since it will be machine-dependent because of the different file manipulation 
procedures employed by different computing sytems. Addition to the data base 
is an easy procedure. 

By way of example, a description is given as Appendix A1 of the retrieval 
system used with the data base at Pratt & Whitney Aircraft. A tape-copy of the 
data base itself has been supplied to the National Aeronautics and Space 
Administration (NASA). 
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Example of Data Bank Input 


CATEGORY: Isothermal Flow SUBCATEGORY: Flow Field Structure 

SOURCE: Habib, M. A. and J. H. Whltelaw: Velocity Characteristics of 

Confined Jets With and Without Swirl. ASME 79-WA/FE-21, ASME Winter 
Annual Meeting, December 1979. (Imperial College) 

EXPERIMENT DESCRIPTORS: 


Type: 

Dimensionality: 

Medium: 

Injectant: 

Description: 


Experiment 

2-D 


Air 

None 

Coaxial jets (with and without annulus jet swirl) 
undergo sudden expansion into cylindrical duct 


1€ASURED PARAf€TERS: 


Mean axial velocity 

RMS axial velocity fluctuation 

VARIED PARAMETERS: 

Annulus to central jet axial velocity ratio 
Annulus flow swirl number 


INITIAL CONDITION MEASURED: 

Axial velocity, tangential velocity, turbulence kinetic energy and 
turbulence momentum transport (LV) 

INSTRUMNTATION: 


Laser anemometer. 
Three-hole pitot probe 

BENCH MARK EXPERIMENT RATING: 6 


COMMENTS: 

Experiments were conducted for three cases: (1) annulus to center jet 
velocity ratio of 3 and zero swirl; (2) velocity ratio of 3 with 
annulus flow swirl number = 0.23 and no center jet swirl; (3) velocity 
ratio of 1.0 and zero swirl. Reynolds numbers of annulus flow were 
77,500, 76,000 and 50,500 for the three cases; Reynolds number of 
center jet was 18,800, 18,540 and 35,500 for the three cases. Seeding 
technique employed was silicone spray separately injected into plenums 
supplying the individual flows. Reasonable agreement between these LOA 
measurements and pitot tube and hot wire measurements is shown. 

A configuration similar to that employed in these experiments was later 
eitployed in experiments conducted by Johnson and Bennett. The Johnson 
and Bennett passage radii differed to some degree and the tube 
separating the inner and outer jet had a tapered trailing edge as 
opposed to a blunt trailing edge employed in these experiments. The 
Johnson and Bennett work includes measurements of (scalar) 
concentrations and concentration fluctuations. Data on the initial 
conditions found in the data bank have been furnished by the authors 
via private communication. 

This work has been published in ASME Journal of Fluids Engineering, 

Vol. 102, (1980), pp . 45-63. Additional information m'y be found in: 

Habib, M. A and J. H. Whitelaw: Velocity Characteristics of a 

Confined Coaxial Jet, ASME Journal of Fluids Engineering, Vol. 101 
(1979), pp. 521-529. 



7.1 Selection 



Benchmark quality test cases were selected from the data base. A total of 13 
test cases were submitted for approval. From this list of cases, eleven were 
selected by the National Aeronautics and Space Administration for calculation; 
(Table 6) several of the selected cases involved multiple operating conditions 
or geometry variations. 

From Table 6 it may be appreciated that Test Cases 1-5 inclusive and Test 
Cases 10 and 11 are flows that are representative of the majority of the 
features contained in the primary zone of gas turbine engine combustion 
chambers; Test Cases 6-8 address features of the dilution zone. Test Case 9 is 
primarily appropriate for afterburners and ramjets. 


Table 6 

Selected Test Cases 


Case No. 

1. Johnson & Bennett -- Coaxial non-swirling suddenly-expanded but 

confined jets 

2. Johnson & Roback -- Coaxial swirling suddenly-expanded but 

confined jets 

3. Habib & Whitelaw -- As ease nos. 1 & 2 

4. Vu & Gouldin -- Confined coaxial swirling streams 

5. Yoon & Lilley -- Confined suddenly-expanded swirling jet with 

outlet restriction 

6. Crabb & Whitelaw -- Single jet in crossflow 

7. Khan & Whitelaw -- Row of jets in crossflow 

8. xnan & Whitelaw — Opposed rows of jets In crossflow 

9. Roquemore et al -- Confined coaxial jets behind bluff body 

{reacting and non- reacting) 

10. Mellor, Chigier & Beer -- Hollow cone liquid spray in air 

11. Semerjian & Segalman -- Radiation (spectral) calculation 



The concentration on nonreacting flows should be noted. In addition, Cases 1-3 
appear to involve excessive redundancy for suddenly expanded Jets, as do Cases 
2-b for swirling flows. This ^Hnphasis arises for three reasons: (1) the 
practical importance of these flows, (2) the deficiencies and limitations of 
the individual experiments, and (3) the importance of thoroughly assessing the 
significance of the turbulence model, the correctness of which is fundamental 
to achieving accuracy in any turbulent flow calculation. 

The reasons for the selection of each case will be discussed briefly. 

7.1.1 Test Cases 1 and 2 - Johnson & Bennett and Johnson & Roback 

These two test cases are really an extension of the same experiment. The 
experiment is shown in Figure 6. Measurements without swirl and with swi**l in 
the outer jet are covered. It satisfies most benchmark criteria, but fails to 
provide known boundary conditions, a tragic omission. This experiment is 
particularly important as a benchmark because it uses a configuration that is 
close to a gas turbine engine (can) combustor, and because of the potential to 
extend the configuration progressively all the way to reacting flows. It is 
also important because of its uniqi.e turbulent mass transport measurements. 

7.1.2 Test Case 3 - Habib & Whitelaw 

The configuration of the experiment, Figure 7, is very close to the 
configuration of Test Cases 1 and 2, and again, both swirling and non-swirling 
cases are covered. It therefore provides an interesting backup to these cases. 
As a benchmark the experiment satisfies some criteria to a degree; it fails to 
provide adequate boundary conditions. Comparison measurements using a pitot 
tube, hot wire probe and LDV (I'iser doppler velocimeter) are included, 

7.1.3 Test Case 4 - Vu & Gojldin 

This case covers similar ground to Cases 1-3, and represents a further 
co-axial jet system; however, both jets are swirling and there is no sudden 
expansion. Of particular interest is the fact that the swirl direction of the 
outer jet was reversed to provide data for both counter swirling and 
co-swirling jets. Although this experiment ha:> been viewed as one with 
progression in flow complexity to reacting flows, this is not really the case 
since a different swirl system was used for tne combusting cases. Measurements 
were made using hot wires. Figure 8 gives the configuration. 

7.1.4 Test Case 5 - Yoon & Lilley 

This is a further swirling flow experiment using a single swirling jet and a 
sudden expansion. It s lulates a gas turbine engine can combustor as indicated 
in Figure 9. A contraction is provided at the outlet to minimize reverse flow 
of external fluid along the centerline. There may also be some difficulty with 
flow separation off the flat swirler vanes used. Measurements were made with a 
5-hole pitot tube. As a benchmark experiment it satisfies some of the 
criteria; however, only one sub-case is to be considered. 
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Figure 6 Johnson et al, Cases 1 and 2 (References 35-36) 
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Figure 9 Yoon & Lilley, Case 5 (Reference 39} 


7.1.5 Test Case 6 - Crabb & Whitelaw 

This is a single jet in crossflow. It satisfies a few benchmark criteria, 
although there are some difficulties with the boundary conditions. It pi'ovides 
an introduction to the multiple jet experiments to follow. 

7.1.6 Test Cases 7 and 8 - Khan & Whitelaw 

These test cases are an extension of Test Case 6, and of each other. They 
involve multiple jets. Figure 10 illustrates the experiment. The data are 
primarily based on pitot tube measurements. 
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Figure 10 Khan & Whitelaw, Cases 7 and 8 (Reference 42) 


7.1.7 Test Case 9 - Roquemore et al 

This experiment is illustrated in Figure .i and involves a reacting flow. The 
flow is not like that in a regular gas turbine combustor. However, it is an 
experiment that is backed by extensive isothermal measurements and flow 
visualization. Unfortunately, it is a flow that is nonstationary. It also has 
very low heat release levels. With respect to the benchmark criteria it is 
otherwise generally a good experiment. 

7.1.8 Test Case 10 - Mellor, Chigier & Beer 

The experiment is far from ideal and is very limited in scope. However, in 
view of the limited availability of such experiments, there was little choice 
but to include it. The liquid spray is of water, and measurements were 
obtained photographically. Figure 12 shows the experimental arrangement. 
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7.1.9 Test Case 11 - Semerjian & Segalman 


APiftINAL PAGE IS 


This case is a numerical test case for direct comparison of the radiation 
model with a calculation using a substantiated and more exact method. 

7.2 Review of Test Cases 

A test case may inply a series of experiments conducted by a particular group 
of researchers at a specific laboratory. 

The most common shortcoming encountered almost universally in the test cases 
was found to be the lack of well-defined inlet boundary conditions. The 
importance of these quantities has not been fully appreciated. 

The computer codes require as inlet boundary conditions the following 
information concerning the entering flow, in addition to obvious property 
values: profiles of mean velocity components, profiles of turbulence 

intensity, and a characteristic dimension Ij, related to turbulent eddy size. 

The turbulence intensities are used to calculate the specific kinetic energy 
of turbulence K, 


K = 1/2 (u'2 + v'2 + w'2) 


(61) 


where u' , v' and w' represent the fluctuating velocities along the x, r,0 
axes respectively. The characteristic dimension is used to calculate the 
dissipation rate of turbulence energy e. 


e = 


Cu(K) 


3/2 


(62) 


where is a constant cormonly taken as 0.09. 

The turbulence model yields the eddy viscosity as 

- 

The Reynolds stresses are needed for closure to the Reynolds-averaged 
equations of motion, and they are then formed from the stress-tensor. 


9u. 

- pu' .u' . = H- (— r — + 
^ 1 J t 9X 


du. 


i i 


for variable density flows. 

(Away from the inlet, transport equations are solved for K and e. ) 


(63) 


( 64 ) 


SO 



If the initial mean velocity components are not adequately specified, the 
equations of motion will result in the calculation of unrepresentative flow 
trajectories immediately downstream of the inlet plane. If the inlet 
turbulence quantities are also incorrectly specified there is an additional 
small but apparent effect on the flow trajectories through terms obtained from 
Equation 64 and containing Mt 3rid K (see Figure 2). 

Incorrect specification of inlet turbulence quantities can have a further 
adverse effect if the flow is a reacting one. This is because the eddy 
break-up combustion model relates the reaction rate of fuel to the eddy 
lifetime, K/c. An incorrect specification hence results in an incorrect 
density field. Since density is strongly coupled to the flow field, an 
incorrect density changes the flow field completely. The incorrect density, 
especially if it is lower than the actual density, cannot be corrected in the 
downstream region as the length scale is modified by the internal flow field. 

The importance of inlet conditions to solution accuracy is highlighted in 
Reference 70. 

A number of pertinent swirling flow cases were selected. Most of these were 
bedeviled by inadequate specification of inlet boundary conditions, and 
several important features of swirling flows were also not satisfactorily 
covered by these experiments. These features were: 

1) Delineation of conditions whereupon a swirling flow is inertia 
dominated as opposed to being turbulence dominated. 

2) Exploration of the critical swirl number range for confined flows, 
including the influence of the confining geometry shape and the 
actual distributions of velocity in exiting the swirl generator. 

3) Effects of backpressure on confined swirling jets. 

The work of Lilley and coworkers (References 39, 71-74) does look at some 
aspects of items 2 and 3. However, the instrumentation used is intrusive and 
yields time-mean information only. Intrusive instrumentation can change the 
nature of the flow being measured, and time-mean information alone is 
insufficient to completely test turbulence models. Furthermore, the number of 
traverse stations used in these experiments is severely limited. This does not 
permit a comprehensive flow picture to be built up. In addition, under certain 
circumstances the flow suffers from formation of a central vortex core 
extending upstream from the exit. The flow is therefore not specified by the 
measured inlet conditions. The swirl generator used has flat vanes, and at 
high angles of attack these vanes are stalled. The hub of the swirler is large 
and has a bluff trailing edge. The generator itself is situated close to the 
dump plane. The swirling flow therefore has a fairly high degree of 
circumferential nonuniformity associated with it, and perhaps some radial 
instability. These conditions are not likely to be reflected in the measured 
inlet conditions. 
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Although the Lilley study is yielding useful qualitative information on 
swirling flows confined in sudden expansions and is highlighting areas for 
further investigation, in current form it is not really suited for the present 
purposes of model validation. 

Yoon and Lilley (Reference 74) have observed that backpressure on a confined 
swirling jet can dramatically influence the nature of the recirculation 
established, Figure 13. It should be noted that the relative dimensions in 
Lil ley's experiment with exit blockage are close to those of a real combustor. 
This flow and effect are therefore clearly of great practical significance. 
Similar effects have been reported in furnaces by Mathur & Maccallum 
(Reference 75) and by Baker et al (Reference 60). This effect needs to be 
explored in detail, particularly in relation with combustor dome shape. 
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Figure 13 Influence of Back-Pressure on Swirling Flow (Yoon & Lilley) 
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It would be extremely useful if isothermal and reacting flow experiments could 
be carried out in the same apparatus so that a natural progression in 
understanding could be followed. Reacting flows alone invoke simultaneous use 
of several of the models and separation of effects can be difficult. 

Generally, such progression has not been done. Gouldin for example, made 
extensive isothermal flow measurements in one apparatus (Reference 38) that 
was chosen as Test Case 4. However, the experiment was extended to reacting 
flows in a different apparatus (References 76 and 77) using a different 
central swirl generator. Although some cold flow measurements are given in the 
new apparatus (Reference 77), they are not extensive enough to be useful. The 
"same" flow field is also being studied by Sommer (Reference 78), but using a 
central swirl generator that is different yet again. Preliminary experimental 
results at comparable swirl numbers are compared against numerical 
calculations made by Ramos for Vu & Gouldins' apparatus and against Vu & 
Gouldins' measurements. The measured results and calculations are compared as 
if the two sets of apparatus were identical, even though the measurements show 
the flows are cT»arly very different, especially for the case of co-rotation. 

When an experiment does progress from isothermal flow to reacting flow the 
progression is usually done in a single step. The intermediate steps of 
isothermal constant density multiple specie, variable density multiple specie, 
or variable density single specie, nonreacting flows are neglected. A reacting 
flow simultaneously exercises the turbulence and combustion models at the 
least. It should be remembered that the turbulence model covers transport of 
both momentum and mass. Since the mass transported could be fuel, it is 
important that this aspect of the turbulence model be explored separately 
before seeking to evaluate the combustion model. 

Few experiments in the literature provide sufficient measurements to permit a 
separate evaluation of the turbulent mass transport in the flowfield prior to 
considering the reacting flow condition. Exceptions to this generality are the 
experiments of Johnson & Bennett (Reference 35), Johnson & Roback (Reference 
36), and Roquemore et al (References 43 and 79). 

Combust ing-f low experiments (as opposed to industrial tests) commonly t'?nd to 
be limited in heat release rate for three reasons: 

1. facility limitations, 

2. operating costs, and 

3. mensuration difficulties 

As a result, such laboratory experiments are frequently carried out at 
atmospheric pressure, modest to low inlet air temperatures, and low overall 
equivalence ratios. The consequence of such operating conditions is that 
combustion is almost always reaction-rate limited. 

The combustion and radiation sub-models chosen for the calculation procedure 
are as simple as possible. As the models were to be simple, the choice was 
made appropriate to the major application, i.e., the modern, efficient, 
high-power gas turbine engine combustor, operating at greater-than-idle speed 
conditions. 



It becomes extremely difficult to find suitable laboratory experiments against 
which to test such models. The actual test cases used were reaction rate 
controlled and had such low radiation fluxes that they were actually "noise" 
compared to the fluxes In the application. 

The test case selected for combustion and radiation model comparisons was Case 
No. 9 in Reference 78. This choice was made because of the Isothermal, 
variable density companion measurements available that permitted some 
amelioration to be made of the difficulties of getting the turbulent mass 
transport correct. Unfortunately, the reacting flow measurements were made at 
atmospheric pressure, ambient inlet air temperature, and equivalence ratios 
less than 0.1 (due to a rich blow-off problem. Reference 52). As was expected 
under such circumstances, comparison of temperature profiles revealed the 
effects of a finite reaction rate that is just not present in the model. 
Similarly, the maximum radiation flux comparisons for this experiment were 
made at a level less than one-twentieth those of a gas turbine combustor. 

The liquid fuel spray model has a number of facets which have to be examined: 
droplet tracking, turbulent diffusion of droplets, droplet heat-up, fractional 
distillation of a multi-component fuel, and coupling of the gas and liquid 
phases. The experiment selected to evaluate the spray model (Test Case No. 10, 
Table 6, and Reference 45) enabled study of only one of these - droplet 
tracking in a uniform airstream. This of course, was a desirable feature for a 
first step in a liquid system in that several simultaneous processes did not 
have to be separated. Water was injected and atomized in a low turbulence, 
ambient temperature airstream at a mass ratio of 4 x 10"3, water to air, and 
at this loading the water droplets did not significantly influence the flow 
structure and did not undergo evaporation of any significance. 

The final constituent flow to be considered was a three-dimensional one, the 
jet-in-crossflow, as both a single jot and as a line of jets. The experimental 
data for this class of flows is generally in fair shape. Difficulties arise 
concerned with representing initially circular jets with the constraints of a 
rectilinear grid, attempting to control numerical diffusion, (a hopeless task 
with present computer storage), and the cost of the calculations. There is 
confusion concerning the operating conditions for the lines of opposed jets, 
(Reference 42), that prevents meaningful comparisons of measurements and 
calculations from being made. Comparisons of concentrations for the line of 
jets cannot be made as the experiment is non- symmetrical. 
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8.1 Considerations for Evaluating Model Performance 
The models studied were for: 

1. Turbulence (momentum and mass transport) 

2. Combustion (for heat release) 

3. Liquid fuel spray 

4. Gaseous radiation 


These physical models are not independent of each other. The interactions 
between models are strong, and this coupling is illustrated in Figure 14. The 
interrelationships between models dictate that a hierarchy of models be 
established, as Figure 15 shows. Note the appearance in Figures 14 and 15 of 
numerics, and its premier position in the hierarchy. Although the subject of 
numerics is not a physical modeling topic, it dominates all the subsequent 
physical modeling through the false diffusion it introduces via the generic 
TEACH solution procedure. Although numerics is not part of the study 
undertaken, its effects must be constantly remembered when considering the 
performance of the physical models. 


Numerics 


Turbulence 


Fuel spray 



Figure 14 Interrelationship of Physical Models 
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Figure 15 Hierarchy of Physical Models 


Figure 1 gives a diagram of the calculation process. The position of the 
physical modeling section is clear and the influence of the following section 
has been described. One further section is relevant - the assembly of the 
equations. 

When the modeled conservation equations were assembled, they were written in 
unweighted form, and terms involving fluctuating density were neglected. The 
implication of this step is that density differences in the flow are not large 
- obviously not too good an assumption for combusting flow cases where 
temperatures in the field can vary widely from point-tc-point and with time. 

For variable density flows a density-weighted decomposition of the 
conservation equations is sometimes used. This is known as Favre averaging. 

For example, the instantaneous velocity in unweighted form is written, 

u = u + u' (65) 

where u and u' are the time-average and fluctuating velocities respectively, 
and in density-weighted form, 

Pu = ^ + pu" = pu + pu" (06) 

where the tilda denotes a Favre-averaged velocity, so 
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The choice cf v^eighted or unweighted equations is somewhat arbitrary 
(Reference 13). The Favre-averaging approach requires modeling of terms like 
p'u' , which then appear in the equations. For these reasons (Reference 12), 
Favre-averaging was not used. Hov/ever, use of Favre-averaging yields 
density-weighted property values. Some instruments measure weighted values 
while others may measure unweighted values. This should be borne in mind when 
considering coi.iparisons of calculations and measurements. 

As will be realized from the cautions presented above, unambiguous assessment 
of the performance of the physical models proved to be very difficult. There 
is strong evidence to believe that many of the shortcomings which in the past 
have been blamed on the turbulence model typically, were in fact, to large 
measure due to poorly specified inlet boundary conditions (Reference 70), and 
to numerical diffusion. 

8.2 Calculation of Test Cases 

The results of the calculations of the test cases will be presented by flow, 

i.e. : 

1. coaxial jets, 

2. swirling jets, 

3. variable dene y Jets, 

4. liquid spray, 

5. reacting flow, 

6. jets-in-crossflow, and 

7. radiation. 

Only a small portion of the total of the calculations and comparisons of 
calculation and data can be made. Results will be chosen to illustrate 
specific points. Additional information was provided to NASA through specific 
task reports. 

8.2.1 Coaxial Jets 

The study of coaxial jets is actually a study of coaxial jets in confined 
sudden expansions, using Test Cases 1 and 3 in the configurations shown in 
Figures 6 and 7. The physical model exercised is the turbulence model. Case 1 
uses water-flow and Case 3 uses air-flow. Table 7 gives the relative 
velocities and dimensions, information on the type of inlet boundary condition 
used, and the grid sizes upon which the calculations were made. Note that the 
relative dimensions of the two experiments were about the same, and that there 
is a progression in jet velocity ratio from 1 to 3. 

Test Case 1 had been predicted blind previously using a combination of 
measured and standard guessed inlet boundary conditions, (Reference 80). The 
predictions had generally been good w'th one notable exception - the 
centerline mean axial velocity variation in the initial region was incorrect. 
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Coaxial Jets in a Confined Sudden Expansion 
(Turbulence Model) 


Experiment 

u /u . 
ann pipe 

*^ann^'^o .pe *^duct^'*ann 

Boundary/Conditions 

Grid 

Habib & 
Whitelaw 

1.0, 3.0 

2.76 

2.81 

Measured 

30x39 

Johnson & 
Bennett 

3.11 

1.93 

2.07 

Measured -i- 2 sets est. 

52x38 


It was felt that the difficulties with centerline axial velocities were due to 
the guessed inlet boundary conditions for the inner jet. An attempt was made 
to "inprove" the assumptions for the inner jet region: 

(i) the dividing lip of the central pipe was modeled as being 
infinitely thin rather than the 0.4 mm used in Reference 80. 

(ii) the turbulence intensity tor the outer jet was improved by 
introducing an increase through the lip-wall boundary layer, as 
shown in Schlichting (Reference 81). 

(iii) the grid was refined to try to capture the potential core of the 
central jet (52 x 38 grid). 

(iv) the diffusing flow at the exit of the central jet was calculated 
using a validated, parabolic matching computer code (Reference 82). 

These conditions were felt to represent improved realism for the input 
boundary conditions over the original assumption, even though any possible 
elliptic effects of the downstream flow on the central jet are neglected. 

Figure 16 compares the measurements of centerline axial velocities with the 
original predictions from Reference 80 and the calculations made using the 
revised boundary conditions for the central jet. Although the present 
calculations do not agree perfectly with the measurements, they apparently 
represent a signif ic:^.nt improvement over the original results. Note that the 
effect extends to at least 10 duct radii downstream from the plane of the 
sudden expansion. 
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ORIGINAL PREDICTIONS 



Centerline Mean Axial Velocities 
Johnson & Bennett's Coaxial Jets in Confined Sudden Expansion 


Figure 16 Revised Boundary Conditions Improve Centerline Velocities 


Unfortunately, the change in boundary conditions destroyed the relatively good 
agreements obtained in Reference 80 for all other quantities at all axial 
stations. An example of this is shown in Figure 17 where a comparison is 
presented of radial profiles of fluctuating axial velocities at 0.836 duct 
radii downstream. The corrected calculations have the jets spreading too 
rapidly compared to the measurements. Too much turbulence is being calculated 
as a result of the new boundary conditions, this is reflected in the turbulent 
transport, and hence in the spreading rates. 

It is interesting to note that when the experiment was originally set up, the 
flow downstream of the sudden expansion plane was not quite symmetric. To 
inprove synmietry, a perforated- pi ate restrictor was introduced upstream in the 
supply pipe for the central jet. Prior to introduction of this restrictor, 
measured centerline axial velocities did fall on the originally predicted 
curve with the original (Reference 80l~5oundary conditions. 
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Original boundary conditions 





Y/Y-MAX 

Alternative boundary conditions 


Example: Fluctuating axial velocity at x/Rq of 0.836 


Figure ]7 Revised Boundary Conditions Worsen Radial Profiles of All Other 
Variables 

The results of the comparison indicate that this experiment is closely 
controlled by subtle inlet boundary conditions. The latest boundary 
assumptions appear less "exact" than the original assumptions. The real 
boundary conditions may fall somewhere between the two sets of assumptions, 
but it is impossible to tell. 

The ability of the code to correctly calculate the effects of jet velocity 
ratio on centerline axial velocities for nonswirling axial jets is illustrated 
for Case 3 in Figure 18. Measured inlet conditions were available for this 
experiment. The quality of radial profiles is illustrated in Figure 19 for 
mean axial velocity at 7.7 cm from the expansion plane for Case 3 at a jet 
velocity ratio of unity. 



■ Calculations 

O • □ ■ Experimental data 

Habib & Whitelaw 

Figure 18 Effect of Jet Velocity Ratio on Centerline Axial Velocity 
Distribution is Correctly Calculated 
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V / 0-0625 

NORMALIZED RADIUS 


Figure 19 Example of Profile Agreement for Test Case 3, (Habib & Whitelaw), 
Nonswirling Jets 


The type, size and strength of the recirculation zones generated are important 
features of the flow. Table 8 compares measured values with those from the 
calculations as a function of velocity ratio. For the velocity ratio of 3.11, 
results for the original (Reference 80) boundary conditions are presented. The 
agreements -.re generally encouraging. 

For Case 1 fluorescein dye was introduced as a tracer into the water flow of 
the central jet to give constant density variable specie mixing of the two 
concentric jets. This permitted the turbulent mass transport side of the 
turbulence model to be investigated. 
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Table a 


Recirculation Zone Comparison for Cases 1 and 3 




Measured 



Calculated 


ann 

u . 
pipe 

Length 

cm. 

Strength 

m/s 

Type 

Length 

cm. 

Strength 

m/s 

Type 

1.0 

36.25 

-5.54 

step 

34 

-4.78 

step 

3.0 

40.0 

-7.56 

step 

34 

-8.04 

step 

3.11 

28 

-0.27 

step 

26 

-0.33 

step* 


♦Original boundary conditions 


The original calculations (Reference 80) supported the reduced value for 
turbulent Schmidt number of 0.5 chosen as a result of an earlier recirculation 
zone study, (Reference 17). This is demonstrated in Figure 20 which cor.ipares 
calculated radial profiles of mean concentratior. of dye for a range of 
turbulent Schmidt numbers with the measurements. However, consideration of 
Figure Z1 in which the calculations are based on the revised inlet boundary 
conditions, indicates the turbulent Schmidt number should be higher than 0.5. 

Case 1: Johnson & 3cnnett 

• UTRC ir DATA 
--- TEACH PREDICTIONS 



Figure 20 Influence of Turbulent Schmidt Number on Mean Mixture Fraction 
Profiles for Original Boundary Conditions 
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Case 1: Johnson & Bennett 


x/Rq ■ 0.836 



ALTERNATIVE BOUNDARY CONDITIONS B 




RADIUS, r/Ro 


Figure 21 Influence of Boundary Conditions on Apparent Turbulent Schmidt 
Number 

Figures 20 and 21 show that opposite conclusions concerning the physical 
modeling can be reached depending on which set of inlet boundary conditions is 
used. 
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8.2,2 Swirling Jets 

The study of swirling jets involved suddenly-expanded and confined coaxial 
jets with swirl in the outer jet (Test Cases 2 and 3), confined coaxial 
swirling streams with both co-rotation and counter-rotation (Test Case 4) and, 
a confined suddenly-expanded swirling jet with outlet restriction (Test Case 
5), Table 9 gives the jet velocity ratios, relative dimensions and swirl 
numbers for these cases; information on the types of boundary conditions used 
and the calculation grids is included. Test Cases 2 and 3 were extensions of 
the nonswirling jet cases discussed above. 


Table 9 

Confined Swirling Jets and Coaxial Swirling Jets in 
Confined Sudden Expansions (Turbulence Model) 


Expertinent 

“ann/'^pipe 

‘^ann^'^pipe 

‘*duct/*^ann , 

S.N. 

Boundary Conditions 

Grid 

Habib & Whitelaw 

3.0 

2.70 

2.81 

0.23 

Estimated 

Various 

Johnson & Roback 

2.31 

1.93 

2.07 

0.47 

Measured • 
estimated 

38 X 33 

Vu & Gouldin 

0.68 

0.67 

3.9 

3.9 

1.0 

1.0 

i • 0.58, 0 • O.S4 
1 • 0.49, 0 • -0.51 

Both measured 
* estimated 

67 X 41 
67 X 41 

Yoon & LiJley* 

0 

1.0 

2.0 

0.53 

Measured 

50 X 32 

•Outlet constricted, 

contraction ratio • 4:1 






For Case 2, Johnson & Roback, the calculations were started at the first 
measurement plane, which was just downstream of the expansion plane. This was 
done because of the uncertainty associated with the inlet boundary conditions 
in the jets themselves. It is unsatisfactory for two reasons; 

1. it gives up on calculating the expansion, and 

2. the input length scales of turbulence (to calculate dissipation) are 
unknown. 

For the present, the central part of the flow was assumed to have originated 
from the center jet and to have a length scale of 3 percent of the pipe 
diameter. The outer part of the flow was assumed to come from the annular jet 
and to have a length scale equal to 3 percent of the annular height of the 
passage. These assumptions neglect both the large dissipation associated with 
the flow expansion, and the fact that the outermost portion of the flow is a 
recirculation zone. The computational representation of the measured velocity 
profiles forming the inlet conditions is given in Figure 22. 
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Figure 22 (continued) 

Shown in Figure 23 is a comparison of measured and calculated mean axial 
velocities along the centerline. It can be seen that the calculations predict 
the presence of the measured central recirculation zone. Its strength appears 
to.be about correct but its length is underestimated by about 15 percent, 
while the position of the forward stagnation point is underestimated by about 
the same amount. 

The general flow development can be described as follows: The development of 
axial velocity is calculated fairly well to 5,08 cm (x/Rq = 0.833). This is 
just about the end of the measured outer recirculation zone. At a station of 
10.2 cm (x/Rq = 1.67) from the expansion, the calculated and measured 
profiles diverge because of discrepancies in the calculated central 
recirculation zone compared to the actual recirculation zone: The calculated 
width of the recirculation is considerably less than the measured width. This 
difference causes the calculated velocity around the recirculation to be much 
lower than the measured velocity because of the reduced blockage. Further 
downstream at 20.3 cm (x/Rq = 3.33), outside the central recirculation, the 
agreement in velocities is improved as the blockage effects are removed. At 
stations further downstream than 30.5 cm (x/Rq = 5.0), a discrepancy in 
velocities appears in the center of the flow, while the calculated and 

measured velocities for the outer region are in good agreement. 
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Roback 



Boundary 
conditions 
mostly measured 


Figure 23 Comparison of Centerline Axial Velocities for Test Case 2 (Johnson 
& Roback) 


The distribution of swirl velocity is in excellent agreement with the 
measurements to a station of 15.2 cm (x/Rq = 2,49) from the expansion. At 
subsequent downstream stations, there is a progressive deterioration in 
agreement of calculated and measured radial profiles. Over most of the radius, 
the measured tangential velocity is approximately constant and only drops off 
to zero on the centerline from a radius of about 2 cm (r/RQ = 0.33). Further 
downstream than 30.5 cm (x/Rq = 5.0) the calculated flow approximates a 
forced vortex while the measured flow approximates a free vortex in the outer 
region, with a forced vortex core. Figure 24. 

There is a large experimental uncertainty in the measurement of radial 
velocity and a large scatter in the data as a result. The calculations of 
radial velocity are qualitatively correct with fair quantitative agreement 
downstream to 2.5 cm (x/Rp = 0.41) from the expansion. The calculations are 
qualitatively correct at b.8 cm (x/Rq = 0.9E), but severe underestimation of 
the maximum values occurs in the outer regions of the flow. At 10.2 cm (x/Rp 
= 1.67) there is very little radial activity and the calculated and measured 
profiles are in good agreement. 
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Figure 24 Test Case 2 (Johnson & Robai'k) at 30.5 cm from Expansion Showing 
Discrepancy Between Calculated and Measured Profiles of Tangential 
Velocity 


The calculations of radial velocity profiles go astray upstream of the axial 
position where the character of the tangential velocity profiles departs from 
that of the measured profiles (between 2.5 to 5.8 cm, and 15.2 cm 
respectively). The discrepancy in tangential velocity profiles is probably 
generated in part by incorrect calculation of the upstream radial equilibrium. 
The measurements of mean axial and radial velocities at downstream stations 
indicate that the tangential velocity profile at 40.5 cm (x/Rq = 6.64) is in 
equilibrium. A free vortex is not a stable profile, and viscosity effects at 
the core modify the profile to the stable forced vortex form locally, as 
indicated at 40.5 cm. Therefore the discrepancy in tangential profiles could 
also be due in part to an inadequate calculation of eddy viscosity across the 
radius. 
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The inlet boundary conditions available for Test Case 3, Habib 4 Whitelaw, 
were less than complete. No details of the swirl generator are provided to 
permit jet interaction effects to be investigated, as was the case for Vu and 
Gouldin's swirling flow experiment, Test Case 4. Inlet boundary conditions of 
distributions of u and w velocities were obtained from Reference 83, as was 
the distribution of kinetic energy of turbulence for the inner passage. No 
information was provided^on K for the outer passage, so some assumptions were 
necessary. In all cases v was taken as zero. 

The first calculation was made on the 30 x 39 grid used for the nonswirling 
case. The turbulence intensity for the outer passage was taken as a uniform 
8.2 percent; length scales for both passages were assumed to be 3 percent of 
the passage height. 

A further calculation was made with the same boundary conditions but the grid 
was coarsened to 20 x 39 to be closer to the grid used by Habib & Whitelaw, 

The results of these two calculations are compared in Figure 25 with the 
calculation of Habib & Whitelaw, and the measured data for centerline axial 
velocities. 
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Figure 25 Centerline Axial Velocities in Test Case 3 (Habib & Whitelaw) 
Compared With Various Calculations 
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Habib & Whitelaw's calculations, made on a 22 x 18 grid, are qualitatively 
correct in that they predict the initial dip in velocity followed by a 
recovery, then a decay followed by a slighi recovery to a downstream value 
that is almost constant. However, the first recovery is underestimated, as is 
the magnitude of the velocity depression that follows it. The downstream 
plateau value is correctly predicted quantitatively. 

The present calculations with the 30 x 39 grid show extremely poor agreement 
with the measured data. Coarsening the grid to 20 x 39 gave an improvement for 
x/De less than 0.5, but the calculations were not even phenomenologically 
correct for x/Dg greater than unity. 

The discrepancies on centerline velocities are so severe that it is not 
worthwhile presenting radial profiles of velocity and other quantities. 

To investigate the flow behavior calculated by Habib & Whitelaw with a 
relatively coarse grid of 18 x 22 gridlines, calculations were made with a 
grid of this density. It was suspected that careful placement of the gridlines 
could be used to selectively supplement with numerical diffusion, the 
diffusion calculated by the turbulence model, as hinted by Jones & Whitelaw 
(Reference 13). 

Unfortunately, the distribution of nodes actually used by Habib & Whitelaw is 
not given in their paper. It was therefore necessary to proceed with their 18 
X 22 mesh in arbitrary fashion. The radial distribution was fixed (and 
arranged to define the passages adequately); the axial distribution was 
initially given an expansion factor of 1.10, then made more uniform by 
reducing the expansion factor to 1.05. Finally, two axial fluid blocks were 
used; each block was given an expansion factor which would produce a dense 
grid at their common face, which was placed in the region of measured maximum 
centerline velocity decay rate. The object of these gridline distributions was 
to examine reduced numerical diffusion in the measured region of rapid 
velocity decay, and allow it to increase in the measured velocity recovery 
region. 

Figure 26 shows the results of this study and the measured data. The results 
are presented in terms of centerline mean axial velocities. The grid with an 
expansion factor of 1.10 has sufficient gridline density close to the 
expansion plane to capture the measured initial depression of velocity. 

However, the increasing numerical diffusion with increasing distance reduces 
the gradient of axial velocity quite severely after about 7 cm downstream. The 
nearly uniform grid line distribution does not capture the initial velocity 
depression; however, the velocity gradient after 7 cm is increased slightly. 

The mesh based on two axial fluid blocks with expansion factors which result 
in a clustering of gridlines in the region from 10 to 40 cm, really increases 
the velocity gradient for distances greater than 10 cm downstream and 
depresses the velocities below those of the other calculations for 30 to 60 cm. 
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Figure 26 Test Case 3 (Habib & Whitelaw) With Coarse Grid Showing Result 
Manipulation Through the Calculation Grid 

Although the study shown in Figure 26 does not reproduce either the measured 
data or Habib & Whitelaw' s calculation, the following can be deduced: (1) 

Habib & Whitelaw' s calculations cannot be grid- independent as claimed; (2) 
placing additional gridlines in the region from 8 to 15 cm with a consequent 
reduction in the region from 15 to 50 cm, might well reproduce the correct 
phenomenological behavior for this experiment; (3) by careful placement of 
gridlines solutions can be manipulated, accidentally or otherwise, to cover up 
shortcomings in the physical turbulence model. 

The experiment of Vu & Gouldin, Test Case 4, provided measured profiles of 
necessary quantities, except turbulence length scale, in the plane of 
confluence for the two swirling streams. However, there was clearly evidence 
in the experiment of an elliptic effect of the downstream flow on the swirling 
central jet prior to it issuing from the central tube, (see Figure 8). The 
implication of this was that it might be better to start the calculation 
upstream of the confluence plane. 
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It was decided that the calculation domain should begin at the downstream face 
of the central swirler for the inner jet, and in the plane of the confluence 
for ^e concentric outer jet. 

The profiles measured at the confluence plane were used for the outer jet, and 
the turbulence length scale was assumed as 3 percent of the jet annular 
height. Inlet boundary conditions for the central jet were calculated. 

The procedure used for the central jet was as follows: The axial swirler is 
fed by means of a long approach pipe of about 19 pipe diameters in length, so 
the flow was assumed to be both turbulent and fully developed at the swirler 
face (stretching the work of Nikuradse, see Reference 81). For the mean 
Reynolds number of the tube the dimensionless velocity profile was calculated 
assuming a power law distribution with an index appropriate to the Reynolds 
number. Since the axial flow swirler has an unswirled flow through its open 
hub, the flow split between the vanes of the swirler and the open port in the 
hub was calculated through matching the one-dimensional flow pressure loss 
characteristics of each separate flow path. The air deviation angles achieved 
by the swirl vanes were estimated as a function of radius, accounting for 
blade pitch to chord ratio, and thickness to pitch ratio. The calculated 
angles varied from 2.9® close to the hub to 63.3° at the tip. From the above 
information the variation of axial and tangential velocities with radius was 
found; radial velocities at the swirler exit plane were taken to be zero, see 
Table 10. 


Table 10 

Central Tube Swirler Exit Conditions 


R cm 

u tn/s 

0.117 

36.82 

0.351 

36.00 

0.539 

35.26 

0.731 

33.17 

0.879 

32.43 

1.039 

31.52 

1.201 

30.39 

1.219 

30.25 

1.321 

29.37 

1.439 

28.12 

1.561 

26.34 

1.653 

24.31 

1.749 

19.75 


e" 

w m/s 

0 

0 

0 

0 

0 

0 

2.90 

1.68 

44.00 

31.32 

59.00 

52.46 

63.20 

60.16 

63.70 

61.21 

65.11 

63.30 

66.06 

63.34 

66 .07 

61.66 

67.17 

57.75 

63.30 

47.40 


ORiGfr'/iL 

HOOR 
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The radial stations given in Table 10 represent the grid lines used to input 
the profiles. The calculated values in Table 10 were scaled as necessary to 
ensure that the step functions actually input to the computer calculation to 
represent the profiles did preserve continuity. 

The assumed values for turbulence intensity and length scale of turbulence 
followed the recommendations of Reference 17 with an intensity of 8.2 percent 
and a scale proportional to the open hub diameter for the central region and 
p'^oportional to the circumferential distance between swirl vanes for the outer 
region. 

The inlet conditions for the central tube were maintained constant for both 
coflow and counterflow swirl cases. 

Figure 27 shows the centerline axial velocity variations with downstream 
distance for counterswirl and coswirl cases respectively. The experimental 
data (References 38, 84) are compared to the present calculations. The zero 
distance position represents the confluence of the two streams. 

Shown also in Figure 27 are calculations of the experiment made by Srinivasan 
and Mongia (Reference 85) and Rhode (Reference 86). Both sets of additional 
calculations were performed with codes that use solution procedures that may 
be described generically as TEACH procedures, differing only in orogramming 
details (Roferences 87, 88, 89). They use the same type of turbulence modeling 
as was used in the present calculations. Therefore, all three sets of 
solutions are directly comparable. The calculations presented should only 
differ from each other through the grids upon which the computations were made 
and the inlet boundary conditions used. 

For grid-independent solutions, which are claimed by all the a. -hors, 
differences exist due only to the boundary conditions. In Reference 85 the 
grid is not quoted directly but is inferred to be 30 x 25; in Reference 86 the 
grid was 30 x 24, while the present grid was 67 x 41. The calculation domains 
were identical for Reference 85 and the present results; for Reference 86 the 
calculations were started at the confluence plane of the streams. 

Consideration of Figure 27 for the counterswirl case shows that all the 
calculations have the same character as the measurements. All utiderestimate 
the strength of the recirculation zone formed and overestimate its length. The 
initial values of the centerline axial velocity differ considerably with the 
two predicted values overestimating the measured value. Surprisingly, the 
initial value for Reference 85 is about twice that of the present calculations 
while the strengths of tne recirculation are comparable. Rhode (Reference 86), 
who used the measurements for the initial values, predicted the weakest 
recirci’l Jition of all the calculations. 

Consideration of the centerline axial velocities for the coswirl condition 
reveals the completely unsatisfactory nature of all three calculations. 


73 


-i , ■ ■■ ‘ 

OF POOR QU^^us I 



S/LU 'AjI00|9A IBWB 9UI|J91U90 



3 C 
O N 
O ’■ 

-, o u 

> 3 ® 

^ O < 

h25 

^ < “■ 

0 (_ Q 

-J < E 
< O < 

1 K O 

g5S; 


“ o 

K S -* C 

= 0. Ul K 

> o ® 
(0 UI O Z 
tt C S “ 
2 O T lU 

I O ^ e 

- S -I = 

t- . < « 

^ w 2 s 
<S S ® 
25 g 

E X 
M E 






O 

« 


o 


o 

<-> 

c 


o 

O 


S/UI 'AI!00|9A IBjXB aU[|J0lU9O 


74 


Figure 27 Comparison of Mean Axial Velocities Along Centerline for Test Case 



OF PCOfV - jy 

The calculations of Reference 85 predict a recirculation that is not only not 
measured but is actually stronger than that calculated for the counterswirl 
case. For Reference 86, no recirculation is predicted but the calculation goes 
completely astray for downstream stations greater than 5 cm. The present 
prediction has the same general character as the other two calculations but 
does not yield a recirculation (just). The initial values of centerline axial 
velocity for References 85 and 86 agree with each other and the measurement; 
the present calculation underestimates this value. 

Figure 27 for ccswirl is not displayed as a performance comparison between the 
three sets of calculations. It does illustrate the differences that can exist 
between three codes of nominally the same solution procedure which use the 
same physical modeling. There is greater difference between the respective 
calculations than there is between any one calculation and the experimental 
data. Such disparities between the calculations must be attributed to 
differences in either grid or assumed boundary conditions, or both. 

A repeat of the present calculations for counterswirl with a coarsened grid of 
31 X 27, which is close to that used by Srinivasan and Mongia (Reference 85) 
for the same calculation domain, increased the initial value of centerline 
axial velocity to about 13 m/s, reduced the maximum reverse velocity of the 
recirculation to about 1 m/s, and moved the forward stagnation point 
downstream and the rear one upstream to reduce the recirculation zone length 
to about 2 cm, see Figure 28. When compared to Table 11 this implies that the 
differences between the present calculations of Figure 27 and those of 
Reference 85 are more due to the inlet boundary conditions than grid 
differences. 
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Figure 28 Centerline Axial Velocities for Test Case ^ Showing Effects of Grid 
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Prediction of Counterswirl Recirculation Characteristics 


Reference 

85 

86 

Present 

Measured 

Forward stagnation, cm 

1.4 

1.1 

0.9 

0.55 

Aft stagnation, cm 

7.2 

6.45 

8.1 

6.0 

Length, cm 

5.8 

5.35 

7.2 

5.45 

Strength, tn/s 

-5.5 

-3.6 

-6.0 

-9.7 

Position max. reverse vel. cm 

2.7 

2.5 

2.7 

3.0 

Velocity grad, from recirculation (m/s)/cm 

1.05 

1.05 

1.57 

2.0 


The actual inlet boundary condition assumptions for the three sets of 
calculations are given in Tables 12, 13 and 14. Table 12 gives the assumptions 
for the outer swirling stream, and compares those for the inner swirling 
stream for Reference 85 against those for the present calculations. In Table 
13 a consistent set of inner swirling stream conditions for the present 
calculations are compared with those of Reference 86. The consistent 
conditions for the present work were obtained by averaging calculated radial 
values at the outlet from the central tube, which was where Rhode (Reference 
86) started f\ .s calculations. 


Table 12 

Test Case 4: Assumed Turbulence Boundary Conditions 


Inner Passage 


Reference 

“‘r^s/u 


It 

cm 

(m^/s^)/s 

85 

4.47 

2.75 

0.0372 

1.1 X 103 

Present 

8.2 

4.2-11.9* 

0.036i-0.088* 

88.5-4.9 X 103’ 

Outer Passage 

Reference 

u'r^c/u 

m^/s^ 

't 

cm 

6 

(m^/s^l/s 

85 

4.47 

1.27-’- 

0.145 

0.0891 X 103+ 

86 

14.1 

14.0+ 

0.0363 

13 X 103+ 

Present 

5.8 

1.3-2. 7* 

0.162 

126.4-252* 


♦Varies as function of radius ■’’Based on nominal conditions 
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Table 13 



Test Case 4: Comparison of Turbulence Conditions at Confluence of Streams 


Inner Passage 
Reference u'nnc/u 


h 

e 


T 

m^/s^ 

cm 

(m2/s^)/s 

Rhode (86) 

14.1 

27.5 

0.0363 

35.9 X 1Q3 

Present 

12.0 

24.6 

0.0313 

35.0 X 103 

(Averaged) 

Assumed Mean 

Table 14 

Velocity Boundary Conditions 


Inner Passage 

Outer Passage 


u 

V w 

1 

V w 

Reference 

ra/s 

m/s m/s 

hi S 

m/s m/s 

85 

30.8 

0 0, 33.18 

calculated* 

0 Calculated* 

86 

Measured 

0 Measured 

Measured 

0 Measured 

Present 

36.8-19.75* 

0 0, 1.7-47.4* Measured 

0 Measured 


♦Varies as function of radius x but close to measured, see Reference 85 


Table 14 contains information on the mean velocities. Note the assumptions of 
zero radial velocity. It is strongly believed that the design of the axial 
swirler would generate significant negative (towards the centerline) radial 
velocity components. 

Table 13 shows that the average turbulence conditions for the inner stream 
came out to be remarkably close for Reference 86 and the present calculations. 
However, the results for coswirling streams are very different, as Figure 27 
shows. This difference is believed to be due to the radial distributions of 
the turbulence quantities actually present in our calculations as opposed to 
the radially uniform values applied in Reference 86. It suggests that a swirl 
number (Reference 90) alone is not sufficient to characterize a swirling flow. 

Figure 29 shows the calculated axial and tangential velocities at 0.2 cm 
downstream from the confluence of the two swirling streams for the 
counterswirl case. Figure 30 displays similar information for the coswirl 
case. Both sets of plots also contain points representing Gouldin's 
experimental data. 
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Comparison of Figures 29 and 30 confirms the assertion that there is an 
influence of the outer stream on the inner stream prior to their actual 
physical confluence at the zero x-station. The central portions of the 
counterswirl and coswirl flows are different in both axial and tangential 
velocities, even though counterswirl was achieved by rotating the outer 
swirler vanes only. 

In both swirl cases the agreement of the predicted tangential velocities with 
the measurements is excellent, except for the maximum values where the 
predictions overestimate the measurements by a considerable amount. The 
agreement of predictions and measurements is fairly good for axial velocities 
for the counterswirl case, although there is some underprediction in the inner 
flow. It is worse for coswirl where the centerline values are in extremely 
poor agreement with the measurements. It is apparent that the calculation is 
going wrong inside the central tube before the two swirling streams interact 
directly. 

There are two possibilities that could account for failure inside the central 
passage: 1) boundary conditions and 2) anisotropic viscosity at the pipe wall. 

It was assumed that the radial velocity component at the swirler exit was 
zero. This is unlikely to be the case. The swirler hub to tip radius ratio and 
the strong variation in air turning angle with radius are likely to induce a 
significant radial flow inward toward the swirler centerline. U would be 
expected that an inward radial velocity component at the swirler exit might 
reduce the excessive outward radial velocity component predicted for the 
central portion of the flow at the 0.2 cm station. This would drive up the 
axial velocities in Figures 29 and 30 and reduce the peak tangential 
velocities to improve the agreement of the initial profiles with the data. 

It is known that a two-component viscosity is really necessary (Reference 91) 
to predict swirling pipe flow correctly. The K-e turbulence model currently 
used does not have this capability. 

Lack of anisotropic viscosity effects in the predictions could represent a 
genuine shortcoming in the model that is adversely influencing the 
relationship of the predictions to the data in the initial plane and the 
subsequent development downstream. However, the shortcomings in the initial 
profiles could equally well be due to incorrect and/or inadequate 
specification of boundary conditions at the exit from the center swirler. 

Figure 31 represents in terms of axial velocities, the downstream development 
of the combined swirling flows for the coswirl case. 
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Figure 31 Downstream Development of Axial Velocity for Coswirling Case of 
Test Case 4 (Vu & Gouldin) 



It can be seen from the experimental data shown in Figure 31 that by an 
x-station of 11 cm asymmetry in the flow has built up, although it seems to be 
reduced again by 24 cm. The predictions and the data agree extremely well at 
axial stations up to 3 cm, in the outer regions of the flow. Agreement along 
the centerline is very poor up to 11 cm from confluence, after which the 
agreement is good out to 24 cm. The outer shear layer is well-predicted out to 
3 cm, while prediction of the inner shear layer although good initially, gets 
worse with increasing distance from injection. 

The equivalent plot to Figure 31 for the counterswirl case shows much better 
agreement of predictions and measurements but has similar weaknesses. 

Although the accuracy of the calculated mean velocity profiles barely 
justifies it, examples of the calculated turbulence quantities are given at 
one specific axial station in the counterswirl flowfield. 

Figure 32 compares predicted and measured radial profiles of axial velocity u, 
tangential velocity w, radial velocity v, turbulence in tens ities u'^n, 3 /u, 
v'rms/iJ» w'rms/i^ J'JSt one of the Reynolds stresses u‘ v7u2, for the 
counterswiri case at an axial station 3 cm from confluence of the two swirling 

streams. Figure 32 gives an idea of the amount of information to be 
assimilated since an equivalent of this figure exists for each axial station 
for each swirl case. 


Th^re was nothing special associated with the choice of the 3 cm station. The 
agreement of predictions with measurements varies with both axial station and 
with the quantity being compared. To present complete profiles of all 
quantities at each station for both swirl cases would result in an inordinate 
quantity of plots. With respect to Figure 32 it should be noted that Gouldin 
has expressed low confidence in the accuracy of the radial velocity 
measurements. 


Also with respect to Figure 32 the following should be noted: the predicted 
fluctuating quantities are derived from the specific kinetic energy of 
turbulence. Equation 4, with the assumption of isotropic turbulence 'uch that, 

K = 3/2 (u'2) (68) 


To calculate the Reynolds stress, the following was used, 



with being obtained from Equation 5. 


(69) 


It can be seen from the figure that the features in the profile of each 
quantity are, in general, well predicted. Numerical accuracy is fair to fairly 
good, but far from being exact. 
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Figure 32 (continued) 
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Figure 32 (continued) 
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Figure 32 (continued) 


To grasp the iniportance of the grid when hybrid differencing is used, the 
present calculations have been repeated on a coarser grid of 31 x 27. Although 
this grid has about the same number of nodes as, for exanple, Rhode's grid, 
the distribution is radically different since the present calculation is 
started inside the central passage. 

The effects of numerical diffusion inside the central passage can be seen in 
Figure 33. Figure 33 compares radial profiles of tangential velocity for the 
two grids 0.2 cm downstream of the confluence of the streams. Smearing of 
peaks and troughs due to numerical diffusion can be identified clearly. It is 
arguable whether Figure 33a) or Figure 33b) represents the better prediction 
in an overall sense . 

In Figure 34 downstream radial profiles of axial velocity for the two grids 
are compared. It can be seen that the coarse grid clearly gives superior 
predictions for both velocity components at the 17 cm station. In this case, 
the diffusion that the turbulence model is clearly failing to provide is being 
supplied by numerical diffusion. 
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b) Hybrid Differencing, 31 x 27 Grid 


Figure 33 Tangential Velocity Profiles Close to Mixing Plane for Counterswirl 
Showing Smearing Effect of Numerical Diffusion in Test Case 4 (Vu & 
Gouldin) 
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b) Hybrid Differencing, 31 x 27 Grid (False Diffusion) 


Figure 34 Demonstration Through Axial Velocity Profiles That Turbulence Model 
Fails to Give Adequate Diffusion in Far Region of Field for last 
Case 4 (Vu & Gouldin) 
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Test Case 5» one experiment of a continuing series by LI Hey and coworkers, 
represented a swirling jet In a confined sudden expansion with an outlet 
nozzle, see Table 9. The swirler angle was selected as representing an 
intermediate swirl Intensity, The nozzle contraction ratio was 4:1 and It was 
positioned such that the duct length to diameter ratio was 2:1. These 
dimensions were considered to be representative of a gas turbine combustor. 

The outlet nozzle also apparently eliminated the core-vortex problem (see 
Section 7.2) for the 45° swirler. 

Measured mean velocity boundary conditions (References 39 and 92) at the plane 
of the sudden expansion were used to start the calculations. The turbulence 
intensity was taken as 12 percent, which LI 1 ley verbally agreed was 
reasonable. The standard length scale assumption of 3 percent of the Inlet 
diameter was used (swirler upstream of expansion plane and adequate geometric 
definition for alternative formulation not being available). The computational 
grid used was 57 x 32 and the calculation domain was extended past the nozzle 
plane In an attempt to ensure the correct back-pressuring effect was 
calculated. Measurements of profiles of mean velocity components were 
available for three axial stations In the simulated combustor. 

In Figure 35 comparisons of measured and calculated mean velocity profiles are 
made for the 15 cm station. It can be seen that the code falls completely to 
calculate this flow, even in qualitative fashion. The failure Is similar at 
the other two stations of 30 and 45 cm. Figure 36 shows the calculated 
stream! ines. 

Figure 13, which is constructed from Reference 74, demonstrates that the 
back-pressuring effect of the outlet nozzle results In the central 
recirculation zone developing a lobed character. The calculations yield a 
conventional recirculation zone shape. Whatever the merits or otherwise of the 
model for swirling flows, it Is apparent that it falls to handle correctly the 
back-pressuring effect of the contraction. 
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Figure 35 Profiles of Mean Velocities at 15 cm for Test Case 5 (Yoon & 
Li 1 ley) 
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Figure 36 Test Case 5 (Yoon & LiTIey), Calculated Streamlines Do Not Show 
the Effects of Back-Pressure on Swirling Flow as Demonstrated By 
the Measurements (compare with Figure l3) 
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8.2.3 Variable Density Jets 

Test Case 9, Roquemore et a1. was used to study variable density jets and to 
exercise the mass transport side of the turbulence model. 

The experimental configuration, given In Figure 11, shows a small central jet 
of carbon dioxide (for present purposes, otherwise propane) with a coaxial 
air- jet that Is removed by the presence of a bluff-body disk flameholder. The 
mixing jets are confined in a long, circular duct. Previous studies of parts 
of this experiment had established that: 

1. the flowfield was dependent on the operating conditions (Reference 
32), 

2. when the central jet was dominant It behaved initially as If 
unconfined (Reference 93), and that a turbulent Schmidt number of 0.7 
seemed appropriate for turbulent mass transfer. 

The study reported in Reference 93 for dominant central jets was directed 
towards operating conditions where the outer air-jet in the experiment had 
almost negligible velocities relative to the CO 2 jet. It was necessary to 
establish If the 0.7 value found suitable for turbulent Schmidt number 
remained acceptable for the central jet dominant flow regime when the two jets 
had more comparable velocities. 

The operating condition was selected as an airflow rate of 2 kg/s and a CO 2 
flow rate of 16 kg/hr. Calculations were made using a range of turbulent 
Schmidt numbers and calculated profiles of mean mass fraction of CO? were 
compared with measured profiles to ascertain the best fit to the data. 

Figure 37 gives the centerline decay of CO 2 for cr^'s of 0.3, 0.5 and 0.7, 
and Figure 38 gives radial profiles of CO 2 at 4 cm from the flameholder for 
this range of Schmidt numbers. The merit of the 0.7 value Is obvious. 

These calculations were made on a 49 x 46 grid, and the grid line placing was 
carefully refined to reduce numerical diffusion to a minimum by ensuring that 
central differencing from the hybrid scheme was used in all sensitive areaf. 

The boundary conditions used for this experiment have been described in 
References 32 and 93. The velocity profile of Reference 93 will be used for 
the unmeasured fuel jet. However, a further refinement has been made for the 
fuel tube turbulence characteristics. The turbulence distribution measured in 
a wind tunnel by Reichardt (Reference 81) formed the basis of these 
characteristics. It was assumed that w'^n^s was equal to v'^ms 
Reichardt 's measurements such that, 

K = 1/2 (i]'2 + 2 v'2) (70) 

and his peak axial turbulence Intensity was scaled to match the centerline 
value measured by Roquemore at his 0.4 cm downstream station. This implies 
that the conditional potential core exists out to this station. These 
assumptions enable the distribution of turbulence specific kinetic energy to 
be determined for the exit plane of the fuel tube. 


92 


PME n 

W 'POOR floRimV 


o 

H 

o 


m 

o 

M 

o 


m 

o 

M 

♦- 

o 



J 






tn 

O 

II 

4-1 

o 



<0 (O 
‘f- w 
■D -M 
<o C 
Oi . 

O I— 

c « 

o t. 

O 

t- 4- 0) 
O) 

■9 E flj 
E 10 L. 
3 (U O 
Z (- E 

4-> 01 
+-> I/) 3 

■a c CT 
■'“5 0 
E o q: 

-C Q ^ 

o 

to E o> 

u 

4 -> (U 

C ^ 1/1 
41 10 

'-■Mo 
3 <0 

-M 

t- C </) 
3 O 01 



4- <0 <+- 


O O 


•M C O 

0 O E 

01 (j -f- 

4- C 0> 
If. O OJ 

UJ O 0 £ 


00 

ro 





J 


94 



ORnMAi Hm n 
OP pocm 


The accuracy of these assumptions for the fuel tube was checked by calculating 
the flow for the conditions at which Roquemore made his downstream 
measurements (1 kg/s air and 4 kg/hr CO 2 ), and comparing predictions against 
measurements at the 0.4 cm station. The results of this comparison indicate 

that the boundary condition assumptions for the fuel tube are indeed good ones. 

It had also been demonstrated in an earlier study of recirculation zones 
(Reference 17) that 0.5 was a more appropriate value for turbulent Schmidt 
nundier than the 0.9 originally recommended, (Reference 10). However, for the 
airjet-dominant flow regime of Roquemore' s experiment use of 0.5 did not 
result in good agreement of the calculations and the measurements. 

Experimental data were available for a case with 2 kg/s airflow and 6 kg/hr 
CO 2 flow, and these conditions represent an airjet dominant flow regime 
where the fuel jet is forced to undergo a stagnation by recirculating air, 

(see Reference 32). Calculations of this case were made with turbulent Schmidt 
numbers of 0.2, 0.3, 0.5 and 0.7. 

Figure 39 ;hows the position of the predicted fuel stagnation point compared 
with the collected data in Reference 43. The prediction does not vary 
significantly with turbulent Schmidt number, and is in perfect agreement with 
the measured data. This agreement gives confidence in this calculation. In 
Figure 40 the effect of turbulent Schmidt number on centerline concentration 
of carbon dioxide is shown, and the results cotnpared with measured data. It 
can be seen that reducing the turbulent Schmidt number improves the agreement 
of the prediction with data. 


It would seem that a turbulent Schmidt number of 0.2 might be appropriate for 
this particular experiment. 
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Figure 39 Fuel Stagnation Points for Test Case 9 (Roquemore et al) 







8.2.4 Liquid Spray 

To evaluate the spray model Test Case 10, (Mel lor et al, Reference 45), was 
selected. The experiment was such that only the droplet tracking aspects of 
the model would be exercised. However, that was considered to be a good 
starting point for a liquid spray which was produced by a real injector. The 
experiment was also valuable in that it provided considerable information on 
the formation of the spray. 

The experiment consists of a 30.5 x 30.5 cm wind tunnel with a convergent 
intake. An air velocity of 14.7 m/s was established, and this was uniform 
across the tunnel. A swir. spray pressure atomizer was mounted on a slender 
sting aligned vertically with th e tu nnel axis to spray downwards. The atomizer 
had a flow number of 5.42 (pph/\/psi) and sprayed 24.57 kg/hr of watei with a 
cone angle of 84.4®. The Sauter Mean Diameter of the spray was measured as 
91.8 See Figure 12. 

To model the experiment the square section wind tunnel was represented as 
being circular in cross-section, with a diameter of 30.5 cm. The atomizer 
dimension was small in comparison with the tunnel dimension so for the first 
part of the spray at least, this difference between the real tunnel and its 
model will be insignificant. The turbulence intensity of the airstream was 
assumed to be 1.5 percent, on the basis of experience with tunnels of similar 
design. The origin of the spray was taken to be the edge of the liquid sheet 
produced by the atomizer. This was given as 0.3 cm downstream from the face of 
the atomizer and at a radius of 0.28 cm. 

The initial boundary conditions for the airstream are rather easily 
established; those for the spray are rather more difficult to obtain. 
Fortunately, sufficient information is provided in the report of the 
experiment to enable good estimates to be obtained for droplet initial 
conditions, if at some effort. 

The optical technique used in the experiment for determinina droplet 
velocities and trajectories enabled the necessary boundary condition 
information to be obtained for droplets at their origin in the ^uel sheet 
breakup. This information was gathered for a large number of droplets, and was 
sorted by drop size. These data were organized statistically to produce 
characteristic velocities and angles as a function of drop size at origin. 

Figure 41 displays a histogram for droplets of 60 microns diameter. The 
distribution of data is nearly normal for the droplet size range 60-90 
microns; however, for both lower and higher drop sizes the distributions are 
skewed, and in some cases suggest two populations. The percentage cumulative 
relative frequencies were therefore plotted in Weibull probability form. 

Figure 42 is the Weibull probability plot for the 60 micron droplets. The 
characteristic velocity (the 63.2 percentile point) in this case is 28 m/s. 
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Figure 41 Histogram for Initial Velocities of 60 Micron Droplets in Test Case 
10 (Mellor et al) 
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Figure 42 Weibull Probability Plot for the Initial Velocity of 60 Micron 
Droplets in Test Case 10 (Mellor et al) 
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The statistical procedure was repeated for each droplet size to produce the 
characteristic velocity plot given In Figure 43. Although there Is soine 
scatter a smooth curve can be drawn through tne data, and this line passes 
through the undisturbed air velocity level for a zero diameter droplet. 



Figure 43 Weibull Characteristic Initial Velocities for Droplets at Breakup 
in Test Case 10 (Mellor et al) 


A similar procedure was adopted for the droplet angle of trajectory at 
breakup. The evidence for there being two populations in the samples appeared 
to be stronger for the case of the angles. However, the accuracy of the 
measurement (taken from the projection of a "double blink" photographic slide 
and based on a line drawn between the centers of two circles only 1-2 circle 
diameters apart) cannot be high, and an increase in the data scatter is 
apparent. In many cases even the Weibull distribution was a poor 
representation of the sample. 

The histogram for 80 micron droplets is given in Figure 44 to illustrate the 
tendency toward two populations in the data sample for angle of trajectory. 

The discontinuity it. the Weibull probability plot of Figure 45 also Indicates 
two populations in the sample. However, for the present purposes the data for 
all drop sizes nave been treated as being a single population. Figure 46 gives 
the resulting characteristic angle of trajectory plot. This plot should be 
treated with some caution due to the uncertainty associated with the data. 
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Figure 44 Histogram for Initial Trajectory Angle of 80 Micron Droplets in 
Test Case 10 (Mellor et al) 
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Figure 45 Weibull Probability Plot for the Initial Trajectory Angle of 80 
Micron Droplets in Test Case 10 (Mellor et al) 
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Figure 46 Weibull Characteristic Angles for Droplets at Breakup in Test Case 
10 (Mel lor et al) 


The effect of turbulent dispersion on the spray is demonstrated in Figures 47 
and 48 for the trajectories of 40, 80 and 200 micron droplets. The 
calculations in Figure 47 have no turbulent dispersion while those in Figure 
48 do. The effect in this case is small, as might be expected foi only 1.5 
percent airstream turbulence intensity. 

Comparisons of predictions with measured data for droplet velocities and 
trajectory angle are given in Figures 49 and 50 at three downstream positions 
in the spray and as a function of drop size. 

The calculations for droplet velocity indicate that by 3 cm from the injector 
droplets with diameters equal to and less than about 35 /«m have reached the 
velocity of the airstream. This indication does not change out to the 9 cm 
position. The measurements indicate that the size for terminal velocity at the 
closest station is about 20 ^m. The agreement of measured and calculated 
velocities is good for droplets of 50 to 200 /'m in diameter at the 3 and 6 cm 
downstream stations. At the 9 cm position the calculations tend to 
overestimate the measured velocities by a small amount. The effect of distance 
is shown for three droplet sizes in the crossplot of Figure 51. The 
calculations agree with the data in shcvHng that droplet velocity increases 
rapidly at all distances with drop size up to about 100 /<m diameter, and then 
increases less rapidly out to 200 jt*m diameter. For distances out to 9 cm from 
the injector a given droplet in the si/e range 70-170 /mn does not lose much of 
its initial velocity. 
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Figure 47 Sample Droplet Trajectories With No Turbulent Dispersion for Test 
Case 10 (Mellor et al) 
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Figure 48 Effect of Wind Turbulence on Di nplet Trajector^’es for Test Case 10 
(Mellor et al) 
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Figure 49 Variation of Average Drop Velocity with Drop Size at 3 Downstream 
Positions in the Spray for Test Case 10 (Mellor et al) 
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Figure 50 Variation of Average Angle of Trajectory With Drop Size at 3 

Downstream Positions in the Spray for Test Case 10 (Mellor et al) 
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Figure 51 Variation of Average Droplet Velocity With Axial Distance From 
Injector in Test Case 10 (Mellor et al) 


The agreements of calculation with measurement for droplet angles of 
trajectory shown in Figure 50 are not as satisfactory as for droplet velocity. 
In general, the calculations overestimate the angle, although agreement seems 
to be improving for droplet sizes above 160 /urn. The measurements show that the 
smaller droplets are, in general, deflected more by the airstream than the 
calculations indicate. However, the 30 ftm drops, which the calculations are 
moving at their terminal velocity by 3 cm, have measured angles in exact 
agreement with the calculated angles at all three axial positions. Despite the 
poor quantitative agreement, the calculations do reproduce the trends shown in 
the data that for droplets up to about 70 /im in diameter the average angle of 
trajectory increases rapidly; for droplets greater than 70 nm the average 
angle tends towards a constant value. 

Figure 52 presents a cause and effect diagram that addresses the difficulties 
evident in Figure 50. The possible causes of the disagreement between the 
calculated angles and those measured are divided into four major categories: 
those concerned with the experiment, the atomization process itself, the 
physical laws incorporated into the model, and the way the data have been 
handled. 
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Figure 52 Cause and Effect Diagram for Figure 50 


Of the causes associated with the experiment, the neglect of swirl and the 
basic techniques are discussed in the paper and the discussion on it. 
Inadequate lighting of the interior of the spray at the liquid sheet breakup 
is felt to be important in that it prevented a complete initial description of 
the spray origin to be obtained. This will be discussed in detail below. 

The droplets in the spray originate from two possible mechanisms: breakup of 
the liquid sheet, and ligament formation from waves formed in the liquid 
sheet. The liquid sheet tends to break up through the development of holes in 

it, and subsequent disruption as sheet continuity is lost. This process tends 

to produce large drops. It would therefore be anticipated that at their origin 
the larger droplets would have velocities and angles close to those of the 
liquid sheet. Figures 7 and 8 of the original paper (Reference 45) and the 
present Figures 43 and 46 support this. The smaller droplets tend to be 
produced by secondary breakup of larger elements generated during the liquid 
sheet disruption, and through breakup of ligaments formed from the surface of 

the sheet by wave motion. The smaller drops hence tend at origin to have a 

wider range of initial velocities and angles than the larger drops. This is 
again supported by Figures 7 and 8 of Reference 45. The liquid sheet having 
two sides, an inside and an outside and often with somewhat different air-side 
boundary conaitions, it could be expected that just downstream of breakup a 
drop sample of small droplets could contain three populations. 


i06 




The nature of this atomization process, together with the inadequate lighting 
of the interior of the spray at breakup precluding measurement of droplets 
originating on the interior with angles less than that of the original sheet, 
would tend to bias the characteristic angles calculated and shown in Figure 
46, to a value greater than the true average for all three populations of the 
smaller droplets. Mel lor et al comment that the small droplets originating 
from the inner liquid sheet surface probably had angles below 20®. Inclusion 
of such a missing population would change the histogram of Figure 44 for 
example, significantly, and would result in a second discontinuity in the 
probability plot of Figure 45 so that three distinct characteristic angles 
could be obtained for droplets of this size. 

Production of three distinct initial angles for a given drop size would 
require a cross-correlation to exist between droplet velocity and angle. Such 
information existed in part (for two of the populations) in the raw data of 
the experiment, but is lost in the data presentation given in the paper. No 
cross-correlation has been attempted in the present work. 

Generation and use of a family of droplet velocity and angle relationships for 
each drop size range as initial conditions for the spray would require some 
changes in how droplet "average" velocity and angles at downstream positions 
in the spray are calculated. Such a procedure is closer to that used 
experimentally and should give better agreement with the measured "averages" 
presented in Figures 49, 50 and 51. Presently, the calculated averaging is 
done through the use of characteristic angles and velocities that do not 
represent all of the droplets originating at the breakup. 

Droplets in the size range of 50-120 /^m appear to have insufficient drag in 
the calculations of their trajectory. Figure 53 shows the drag law used and 
compares it with the data of Yuen and Chen (Reference 94), with which it is in 
general agreement up to a droplet Reynolds number of about 100. Under present 
conditions the critical Reynolds number of 100 applies to droplets of about 90 
fim. Therefore the "dip" in the drag curve would apply for droplets of 50-120 
fm. However, there is no justification to accept this as the explanation for 
the discrepancies of Figure 50. 

Photographs of the spray show that the majority of the droplets are indeed 
spherical . The turbulence of the stream into which the spray is introduced is 
low initially so that turbulent dispersion of droplets does not have a 
significant effect on turbulence, as Figures 47 and 48 illustrate. Neglect of 
modification of the air turbulence by the spray is therefore not significant. 
The spray is a "thin" spray so that drop-drop interactions are believed not to 
be significant. 

The most promising candidate from Figure 52 for the cause of the angle 
discrepancy therefore emerges as the treatment of droplet initial conditions. 
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Figure 53 Drag Coefficient Information for Droplets 


In Figure 54 the variations in mean drop size based on spatial distribution 
with axial and radial position are compared. In this comparison only the 
shapes of the distributions should be compared since different techniques were 
used to obtain these distributions. For the experiment the distributions are 
based on the arithmetic mean drop size obtained at a given distance from the 
real droplets. Because of the wide range of initial conditions with respect to 
angle of trajectory and velocity (see Figures 7 and 8 of Reference 45) the 
droplets are dispersed across the spray close to the injector. In the 
calculations the droplets are assigned a single initial angle and velocity 
according to their size from Figures 43 and 46. The level of mean droplet 
diameter in the calculations was affected by the range of droplet sizes 
considered. This range was 20-180 ju.m and it was divided into / intervals. The 
range was restricted by programming limitations in the computer code for the 
spray model. The actual size distribution ranged from 10-220 /u^m droplet 
diameter. 
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Figure 54 Variation of Mean Drop Size (Based on Spatial Distribution) With 
Axial and Radial Position for Test Case 10 (Mellor et al) 
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It can be seen that the shape of the mean size distribution in the spray at 9 
CPI is well represented by the calculations. This downstream distance 
represents an almost "fully developed" region of the spray where initial 
conditions are less influential and the concept and use of characteristic 
velocities and angles of trajectory for initial conditions in a calculation 
should be becoming appropriate. Thus, at distances downstream from the 
injector the well known migration of smaller droplets toward the center of the 
real spray results in the tendency toward the uniform distribution of mean 
drop sizes shown in the calculations. 

8.2.5 Reacting Flow 

Test Case 9 (Roquemore et al, Reference 43) was selected as the reacting flow 
experiment to evaluate the combustion model. This choice was made because the 
experiment had sufficient depth for the turbulent transport to be essentially 
eliminated as a source of uncertainty in the calculations, see Section 8.2.3. 

The experiment consists of a long axisyrmetric combustion tunnel employing a 
disk flameholder having a cylindrical forebody with a streamline nose. Gaseous 
propane is introduced as fuel from a small orifice in the center of the disk 
while air enters as an annular jet around the flameholder. Figure 11 gives the 
layout of the combustor. 

Unfortunately, the experiment has a serious problem in that it is a 
nonstationary flow, (Reference 53). The nonstationary flow is associated with 
vortex shedding off the bluff body flameholder, (Reference 95), and in the 
range of flow rates of interest the eddy shedding could be acoustically 
coupled with the tunnel natural frequencies, (Reference 96). The nonstationary 
behavior has to be borne in mind when using this experiment. 

The boundary conditions and the calculation grid used were defined in Section 
8.2.3. Since the flow field for this experiment has two distinct (and steady) 
regimes depending on operating conditions (Reference 32), two sets of flow 
rates were used, 2 kg/s air and 4 kg/hr propane (A2P4) for the air jet-dominant 
regime and, 2 kg/s air and 16 kg/hr propane (A2P1C) for the central fuel jet 
dominant regime; an attempt was made to find the breakthrough point where the 
fuel jet broke through the air recirculation, using fuel flow rates of 8-12 
kg/hr propane. For the nonbroken- through regime the turbulent Schmidt number 
was taken as 0.2 and for the broken- through regime a value of 0.7 was used, 
these being based on the variable density isothermal study reported in Section 
8.2.3. The calculations were started in the plane of the flameholder face. 

The calculated flow fields for nonbroken-through and broken-through conditions 
are shown in Figures 55 and 56 respectively. The streamlines of Figure 55a) 
show the multiple stagnation points in the complex recirculation zone 
established in the wake of the flameholder. The main burning zone is attached 
to the flameholder, as Figure 55b) shows. When the fuel jet breaks through the 
air recirculation the main burning zone is downstream, with a pilot region on 
the flameholder. Figure 56b). The simpler recirculation zone shown in Figure 
56a) seems to go with the reduced mixing implied by the higher value of 
turbulent Schmidt number. 
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Figure 55 Reacting Flow Calculations for Nonbroken-Through Flow, Test Case 9 
{Roquemore et al) 
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Figure 56 Reacting Flow Calculations for Broken-Through Flow, Test Case 9 
(Roquemore et al) 
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To ascertain the Importance of radiation, both gaseous and particulate, to tfi< 
calculations the A2P8 reacting flow case was run on the fine grid with a 
turbulent Schmidt number of 0.5 for the following conditions: no radiation, 
gaseous radiation, a “small soot" case and, a "large soot" case. Sooting 
flames were considered because photographs of the flame (Reference 52) 
suggested some luminosity adjacent to the flameholder, and Roquemore 
(Reference 79) mentioned soot depositions on the face of the flameholder 
during the course of the run. 

The "small soot"/" large soot" definitions were as given below. Small soot: 
particle diameter of 250A and number density of 15 x 10° particles per cc; 
Large soot: particle diameter of SOOA and a number density of 20 x 10° 
particles per cc. 

The large soot definition was based on a typical smoke number density/size 
relation derived from the measurements of Stockham and Betz (Reference 97) by 
scaling to a JT9D-7A size, reported hydrocarbon smoke sizes by Whelan 
(Reference 98), and carbon particle size measurements at the combustor exit In 
a JT9D experimental engine. This was felt to represent a likely upper limit. 
The small soot definition was derived from the large soot definition by 
tempering It according to the laboratory flames reported by Magnussen and 
Hjertager (Reference 22). The small soot case was felt to be closer to the 
flame under consideration. Soot, in both cases, was assumed to exist in a 
single zone in the flame to a downstream distance of 9 cm from the flameholder 
and out to a radius of 6.92 cm. The zone was sized based on the photographic 
evidence (Reference 52). 

To make the radiation calculations, the measured temperature distribution on 
the flameholder face was averaged to a value cf 750 degrees K; tunnel wall 
temperatures were assumed to be a uniform 340 degrees K. The emissivities of 
the flameholder and tunnel wall were 0.8 and 0.65 respectively. The former 
value was based on the soot depositions noted for the flameholder and the 
latter was an integrated mean for the tunnel walls and windows. 

Plots of the centerline temperature distributions revealed a maximum 
difference between the no radiation and large soot cases of about 40 degrees 
K, occurring abcut 7.6 cm from the flameholder. Radial profiles of temperature 
at stations 6.3, 7.5 and 25.3 cm from the flameholders gave respective 
differences of 53, 56 and 3 degrees K. The peak difference occurred at a 
radius of 1.3 cm. These temperature differences are very small, and Indicate 
that radiation is not an important source of heat loss in the test case. The 
maximum difference in temperature between the large soot and gas-only 
radiation was about 20 degrees K. This indicates that soot is not important 
either. 

On the basis of the diagnostic testing described above, it was decided that 
gaseous radiation only need be considered. 
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Figure 57 presents information on the position of the fuel stagnation point in 
reacting flow for the A2P8 conditions. The original calculations referred to 
in the figure were made in 1980 using the 39 x 39 grid of Reference 32, o’.!; 
of 0.5 and the original fuel jet inlet boundary conditions. 


original calculations 



Figure 57 Fuel Stagnation Point in Reacting Flow for A2P8 in Test Case 9 
(Roquemore et al) 


The effects of the refined grid and revised fuel jet boundary conditions can 
be seen to be small. The reduction of turbulent Schmidt number to 0.2 moves 
the stagnation point significantly downstream, thereby in^)roving the agreement 
with the measurements. However, breakthrough at these conditions is still not 
obtained from the calculation, whereas the measurements show the flow is 
actually broken-through. 

The flow behavior is such that A2P8 is just about at the combusting case 
breakthrough point. It is a little optimistic to expect to be able to 
calculate such a singularity exactly. Even in isothermal flow the calculations 
required more fuel flow for breakthrough than the experiment indicated 
(Reference 32). 

The earlier calculations had indicated that breakthrough took place off-axis 
due to distortion of the air recirculation zone by off-axis burning (Reference 
96). Examination of the streamlines predicted by the current calculations 
suggested that sufficient off-axis burning was not taking place to give such a 
breakthrough, so that although the fuel stagnation point agreed with previous 
results, the character of the predictions had actually changed slightly from 
that of the original predictions. The change in character could be attributed 
to the reduction in diffusion due to the increased grid density. 
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To be sure that this was indeed so, a coarser grid of 39 x 39 was run with the 
same fuel jet boundary conditions. This was compared with the results of a 49 
X 46 grid calculation; both calculations had no radiation so that the effects 
of grid alone could be studied. It was confirmed by comparing the two sets of 
calculated isotherms. The reduced false diffusion (of fuel) by the finer grid 
moves the peak temperatures of the flame radially inward and upstream compared 
to the coarse grid. 

The unfortunate dynamic character of the experiment must be considered also 
with respect to the breakthrough point. 

The nonstationary nature of the flame might have resulted in the necked down 
air recirculation zone (see Figure 55b) being fractured earlier than would 
have occurred otherwise (Reference 96). The fragile character of the hot gas 
connection between the pilot flame on the flameholder and the downstream flame 
can be appreciated. Early rupture of this connector could be induced by a 
superimposed oscillating pressure field as postulated in Reference 96. It 
might, therefore, be impossible for the present steady state codes to 
correctly predict the breakthrough point fer this experiment. 

Using the indicated turbulent Schmidt number for the broken-through flow 
regime, the A2P16 combusting case was calculated. 

Figure 58 compares calculated and measured radial profiles of temperature at 
several axial stations. It is clear that the character of the flame is 
correctly predicted as burning off the centerline of the combustor. The radial 
positions of the maximum temperature are predicted fairlv well, as are the 
widths of the flame. Maximum temperatures are considerabV overpredicted 
however. 

The discrepancies seen in Figure 58 between predicted and measured temperature 
levels are up to 700 degrees K. This level of discrepancy exceeds that 
normally anticipated for the physical modeling used. 

There are two important sources of uncertainty; that associated with the 
combustion model, which is 200-300 degrees K (Reference 22), and that 
associated with neglect of fluctuating density, which is about 300 degrees K 
(Reference 99). These uncertainties are comparable. Since the study in 
Reference 22 was carried out neglecting fluctuating density effects, the 
uncertainty in the present case would not be expected to exceed 3C0 degrees K. 
It is apparent that some other effect must be accounted for. 
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Figure 58 (continued) 


In this particular experiment there is a source of uncertainty that can be 
associated with the temperature measurement itself. The overall equivalence 
ratio ^s extremely low, suggesting little carbon should be formed in the 
flame. At these flow conditions for distances from the flameholder greater 
than 8 cm (x/df > 0.57) the measured temperatures are less than 1200°K, see 
Figure 58. Very little luminosity is emitted from non-sooting flames by gases 
at temperatures of 1000°K and less. However, direct photographs of the flame 
have been taken. Reference 52. This implies the presence of temperatures 
higher than those actually measured and/or the presence of free carbon in the 
flame. 
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The flame has been shown to be nonstationary (Reference 53), with luminous 
ball-flames and alternating dark spaces in evidence to distances greater than 
x/df of 5. The presence of carbon particles, which give rise to the 
luminosity of the ball-flames when the overall equivalence ratio of the flow 
is so low, indicates that the fuel and air mixture in the balls, and hence the 
flame temperatures, are much higher than the average. The tei^jeratures 
measured by the thermocouple with its relatively long response time (compared 
to the ball-flame frequency) are clearly average values for the ball-flames 
and the dark spaces between them. 

The development of the ball-flames is described in References 95 and 96. A 
feel for what temperatures the ball-flames might represent can be obtained in 
crude fashion as given below for the furthest station from the flameholders 
for which data are available in the A2P16 case, i.e. x of 16 cm (x/df = 

1.15). The 500 frame per second film clips given in References 52 and 95 
suggest that the ball-flames are actually shed from the recirculation area at 
an x/df of about unity, and the calculations of wake development of Reference 
96 are in reasonable agreement with this observation. 

Using spectrophotometers, Roquemore et al (Reference 53) acquired 
time-resolved CH emission intensity records at a number of axial stations for 
operating conditions of 2 kg/s airflow and 18 kg/hr propane flow. At a 
downstream station of x/df equal to 2.86, where it was easy to distinguish 
distinct luminosity peaks on the traces presented in the reference, the flame 
was found to be present about 74 percent of the time. In Reference 100, for 2 
kg/s airflow and 16 kg/hr propane flow, at this station a typical luminosity 
time history plot indicates the flame is present only about 40 percent of the 
time. While the percentage time that a flame exists at a location depends on 
both the axial location and the operating conditions, the values of 40 and 74 
percent appear to represent reasonable possible lower and upper limits for the 
region of interest. 

Through the use of conditional sampling (Reference 100), it was found that the 
most probable centerline mean axial velocity in nonluminous regions of the 
flow was 31.6 m/s at the stated operating conditions, while that in the 
luminous regions was 42 m/s. The value in the nonluminous regions was equal to 
that in the isothermal flow case at those operating conditions. 

The information thus gleaned from References 53 and 100 can be used to obtain 
a crude estimate of the ball-flame temperature by weighting for the 
nonstationary nature of the flow. Thus, 

T = [TfUfTf + (1-T^) Uiso^isol/u 1 end error. (71) 


Or hUuK i 
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= averaged temperature read by long time constant 
thermocouples 

= percent time a ball-flame is present 

= mean velocity of ball-flame 

= unknown ball-flame mean temperature 

= mean velocity of nonluminous region 

= mean temperature of nonluminous region 

u = mean velocity of reacting flow obtained by unconditional 

sampling 

In the present case u was 39 m/s, while T was assumed equal to 293'"K, the 
nonreacting flow condition. 


The end error is used to match the radial profiles to the temperature outside 
the flame region. It arises because the weighting is not exact. 

The expression for T was solved for Tr using the two values for the 
measured velocities, and the assumed Tiso* corrected profiles were 
matched to a temperature of 293"K in the outer regions of the flow. 


Figure 59 compares the measured data and the calculations for the axial 
station 16 cm from the flameholder with a possible band representing where the 
ball-flame temperatures might fall, depending on the actual intermittency of 
the flame. The figure indicates that if the flame were not intermittent, then 
reasonable agreement with predictions might be expected, subject to the 
uncertainties of the combustion model or the type of weighting used for the 
Reynolds equations. 

Recent CARS measurements by Roquemore (Reference 101) indicate the ball-flames 
have temperatures of around 1900''K, while the dark spaces between have 
temperatures around 400''K. These measured temperatures agree well with the 
corrected data band shown in Figure 59. 

The A2P4 case is far enough removed from the actual breakthrough point that a 
fair comparison of the calculations and measurements may be made. Figure 60 
compares calculated and measured mean axial velocities along the combustor 
centerline. It shows that the right strength is calculated for the 
recirculat ion rone, but that the position of the fuel stagtiation point is 
underestimated and that of the aft (air) stagnation point is overestimated. 
This is entirely consistent with the findings for isothermal conditions 
(Reference 30. (The measured velocities are for nonconditional sampling, but 
this appears appropriate within the recirculation zone.) 
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Figure 60 Centerline Axial Velocities for Test Case 9 (Roquemore et al) at 
A2P4 Conditions 


The radial temperature profile comparisons at a number of axial stations for 
the A2P4 case are given in Figure 61. The respective levels of temperature 
should be viewed with reference to Figure 59, and only the characters of the 
calculated and measured flames considered at this point of understanding in 
the experiment. 

It can be seen from Figure 61 that at all axial stations the overall width of 
the flame is correctly predicted, as is its character of off-axis burning 
initially. However, the calculated flame is burning much more rapidly than the 
measured flame. For example, at 6 cm downstream from the flameholder the peak 
temperature is calculated to be already on the centerline and maximum heat 
release is over, whereas the measurements show the flame still to be burning 
off-axis with higher temperatures than calculated. The radial width of the 
calculated major reaction zone is much narrower and the peak temperatures 
exist at smaller radii than the measurements show, at axial stations to 4 cm. 
At the 8 through 18 cm stations the measured temperatures are all greater than 
those calculated, showing that the actual reaction at these stations is less 
complete. 
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Figure 61 (continued) 
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8.2.6 Jets- In-Crossflow 

This important class of flows was depicted by a series of experiments 
represented ly Test Cases 6, 7 and 8, beginning with a single jet and 
progressively extending to multiple jet systems. Test Case 6 is a single jet 
in crossflow using the experiment of Crabb & Whitelaw; the experiments of Khan 
& Whitelaw are used for a row of jets in crossflow (Test Case 7} and for 
opposed rows of jets in crossflow (Test Case 8). The test configuration is 
essentially the same in all cases with just variation in the jet inlet 
ducting, and is illustrated for opposed jets in Figure 10. 

Thu model exercised is the turbulence model, and the computer code is 
3D-TEACH. In addition to the physical model performance in these flows, the 
effects of representation of circular jets in rectilinear form, and the 
behavior of the hybrid differencing in three-dimensions will be encountered. 

To save on grid lines, and hence on computer storage and calculation time, 
flow symmetry conditions were used wherever feasible. 

For the single jet in crossflow, with a jet to cross-stream velocity ratio of 
2.3, the calculations were performed for one-half of the jet and wind tunnel. 
The calculation domain extended 5 jet diameters upstream of the jet and 24 
diameters downstream, for a total axial length of 30 diameters (76.2 cm); ;he 
half-width of tlie tunnel was 8.5 jet diameters (21.6 cm), while its height was 
11.8 diameters (30 cm). The jet had a diameter of 2.54 cm. 

A 35 x 18 X 16 grid was used. Care was taken to obtain a suitable match of the 
grid lines representing the orifice supplying the jet and the grid lines in 
the wind tunnel. The orifice was defined by 5 axial and 10 transverse grid 
lines. Originally, these were distributed uniformly. Although the grid lines 
themselves were uniformly distributed in the complete grid the relative 
displacement of vector and scalar fluid cells resulted in a nonuniform 
distribution of cells. This caused a convergence problem in the calculation. 
The number of grid lines within the jet initial plane was maintained, but 
their distribution was adjusted to give a uniform distributioti of cells within 
the calculation domain. The calculation grid is shown in Figure 62. 

The boundary conditions were a combination of measured inlet conditions and 
assumptions (Reference 102). The wind tunnel velocity was 12 m/s with a 
turbulence intensity of 0.55 percent, both assumed uniform. The length scale 
was taken as 3 percent of the tunnel hydraulic diameter. For the entering jet, 
the length scale was taken as 3 percent of the real jet hydraulic diameter. 
Turbulence intensity was taken to be the same as that assumed by Crabb, i.e.. 


u' . 

3/2 {—^ 21 ) u ^ = 0.02 


A iiieasured v-velocity profile was available for the jet, and this was applied 
to the modeled jet symmetrically. 
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Figure 62 (continued) 


The representation of the jet origin by a symmetrical stairstep pattern 
results in a jet cross sectional area that is 3.34 percent smaller than the 
actual area, and a modeled jet circumference that is 21. b percent greater than 
the actual circumference. The hydraulic diameter of the modeled jet is 
therefore only 80 percent (2.02 cm) of the real jet (2.B4 cm). 

The calculations were made with the synmetrical inlet velocity profile for the 
jet measured by Crabb with zero crosswind. With the stairstep representation 
of the cross section and the fluid block representation of the inlet profiles, 
the modeled jet has a momentum ratio that is 3.6 percent less than the real 
jet; the modeled mass flow rate is 4.5 percent less. The initial blockage due 
to the modeled jet is that of the real jet (width equal to jet diameter). 

The combination of the reduced momentum ratio for the modeled jet and reduced 
effective diameter, together with the increased initial surface ,area, would be 
expected to result in a reduced penetration compared to the real jet. The 
extent to which the jet mixes with the crosswind would be expected to increase 
as a result of the increased surface area. The numerical diffusion that 
results from the grid limitations is additive to these physical effects. 
Therefore, the calculations could be expected to be seriously in error when 
compared to measurements, irrespective of any deficiencies associated with the 
physical modeling embodied in the code. 

The mass deficiency of the modeled jet compared to the real jet was corrected 
using the variable detisity capability of 3l)-TKACH so that jet initial derisity 
was increased to achieve the true mass flow rate. The jet temperature was 
lowered accordingly. With this change the ratio of modeled to actual jet mass 
flow rates was 0.993, which was considered adequate. The effect of the 
artificial temperature change on jet penetration, development and mixing was 
considered to Be negligible. 

To better understand the development of the jet flow visual izat ion is used. 

The visualization method used is streaklines: the simulation is based on the 
water-tunnel visualization technique, where the illuminateii plane has finite 
thickness, although it is not large. It must be remembered, therefore, that 
the streakline plots in the figures are proj^ections of a three-dimensional 
track onto a two-dimensional plane. 

The x-station referred to is referenced to the beginning of the calculation 
domain: To compare with the calculations of Reference 102 for example, 
subtract 0.1397 m from the present x-values to lino up with the jet origin 
centerl itie. 

Figure 63 shows the axial development on the centerline of the jet. It ca;i be 
seen that the flow is highly three-d imens iona 1 behind the jet at all axial 
stations downstream from the hole edge. Figure 64 , at the (side) edge of the 
jet (z - 0.01 m) shows that the mainstream flows around the jet into the wake 
region. This crossflow is evident at a 7 0 of 1, i.e., half a jet diameter 
beyond the edge of the hole (z = 0.02bm) and at a height of half a jet 
diameter off the tunnel floor, ^igure 6b. 
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Figure 63 Flow Visualization of Single Jet in Crossflow, Test Case 6, (Crabb 
& Whitelaw) Using Streaklines - Axial Development on Jet Centerline 
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Figure 65 Single Jet in Crossflow - Flow Into Wake Region Away From Jet Origin 
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The cross-sectional development of the jet is illustrated in Figures 66, 67 
and 68. 

Figure 66 is for a station at x of 0.1m, or about 1.5 jet initial diameters 
upstream of the jet origin centerline; see Figure 63 for orientation. Two 
important features are apparent: circulation is developing upstream of the jet 
and, evidence of a bow wave can be seen developing. 

Figure 67 illustrates the vortex development in the jet cross-section at a 
station of 0.2m. In this case, only a single vortex pair is visible 
(calculation is symmetrical about the jet centerlineK It is not known whether 
the grid is inadequate to resolve the second pair of vortices or if they are 
physically absorbed into the main vortex pair in this case. The low initial 
jet penetration suggests that the latter explanation might be so, although a 
(rather unrealistic looking) additional vortex pair was calculated by Claus 
(Reference 102) for this experiment using a finer grid. If the calculations 
for this station are compared with those in Reference 102, the closest common 
station is that given for Figure 3c of the reference (x of 0.1746 m in the 
present coordinates). The comparison would be of streaklines in Figure 67 to 
velocity vectors in the reference. 

Jet penetration and spreading can be seen by comparing Figure 67 to Figure 68 
which is at a station of 0.5 m from the origin of the calculation. The centers 
of the vortex pair have moved upwards and transversely outwards. 

F'jure 69 presents profiles of mean axial velocity, normal to the wind tunnel 
flow and in the plane of the jet centerline. The profiles are at a number of 
streamwise positions. Velocity is nondimensional ized by the undisturbed 
cross-stream velocity, and distance by jet initial diameter. 

At an x/D of negative unity, the point of inflection in the calculated 
profiles is noteworthy. It could represent the influence of the horseshoe 
vortex referred to above. The measured approach profile appears to oe fairly 
well-represented by the calculations. On the jet vertical centerline (x/D 
equal to zero), there appears to be disagreement between the calculations and 
the measurements. If real, this could be due to the distortion of the jet 
velocity profile inside the supply pipe by the crosswind. This condition was 
not featured in these calculations due to lack of representative boundary 
conditions. 

Immediately downstream of the jet origin things are changing rapidly. At an 
x/D of unity the calculations show a recirculation at the wall, while the 
measurements show a positive velocity adjacent to the wall. A positive 
velocity immediately behind the jet could be generated through crossflow from 
the horseshoe vortex. If so, the grid used in the calculations is not fine 
enough to resolve this feature. At subsequent downstream stations, x/D from 
2.0 to 14, the wake region appears to be effectively captured by the 
calculations. However, the rapid mixing postulated for the calculations would 
result in rapid dissipation of the excess axial velocity of the jet over that 
of the crosswind. The calculations do yield this result, whereas the excess 
velocity persists in the measurements to at least x/D of 14. 
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Figure 67 Single Jet in Crossflow - Development of Double Vortex Just 

Downstream of Jet Origin (Calculation Sytiretrical About Centerline) 





Figure 68 Single Jet in Crossflow - Jet Penetration and Transverse Spreading 
Through Movement of Vortex Centers 
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Figure 69 (continued) 


Loss of the excess velocity makes it impossible to discern the jet trajectory 
from the calculations. Claus (Reference 102) required a 40 x 30 x 20 grid to 
capture the excess velocity. 

Figure 70 shows isovels at a station x/D of 8 and compares the calculations 
against the measurements taken from Reference 41. The calculated isovels for 
y/D greater than 6 are speculative since the grid is too sparse in this area 
for good definition. In Figure 71 the present calculations are compared with 
those of Crabb. 
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Figure 70 Single Jet in Crossflow - Comparison of Calculated and Measured 
Isovels of Mean Axial Velocity 
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Figure 71 Single Jet in Crossflow - Comparison of Crabb and Present Isovel 
Calculations of Mean Axial Velocity in Y - Z Plane at x/D of 8 
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The nature of the present calculations is in close agreement with Crabb's 
calculations. Crabb used a 22 x 12 x 12 grid with the jet origin represented 
by 6 X 6 nodes. Consideration of Figure 68 helps in understanding Figure 70; 
some direction of flow arrows from streaklines have been sketched on Figure 70 
from the earlier figure. The calculated vortex pair (symmetrical about the jet 
centerline) consists of a central region of low axial velocity with a pair of 
"wings" of high axial velocity wrapped around two-thirds of its circumference. 

Figure 70 shows that the calculated jet development is much advanced over the 
measured development. The central region of low axial velocity is much more 
diffuse radially and the wrap-around of the high axial velocity "wings" is 
greater, with the 1.04 u/U «5 isovel divided into two separate islands on either 
side of the centerline. These results are consistent with the expectations 
from the present modeling of the jet and the limited calculation grid used. It 
is most unlikely that a grid independent solution was obtained. 

Shown in Figure 72 is a cross-sectional comparison at an x/D equal to 8 of 
measured and calculated concentration contours of a helium trace gas 
introduced into the jet. The mass fraction calculations were made with a 
turbulent Schmidt number of 0.5. 

The calculated jet can be seen to be penetrating into the crossflow much more 
than the measured jet. The vertical spreading of the jet is about the same as 
measured (despite the increased penetration); the transverse spreading 
however, is somewhat less than measured, although the calculated 0.2 mass 
fraction contour has already split into two separate closed lobes. This is 
consistent with the jet development presented in Figure 70. 

The present calculation of concentrations gives results that are comparable to 
the computations of Claus (Reference 102) using a 20 x 20 x 12 grid. In the 
plot of Reference 102 the calculated contour for 0.2 mass fraction of helium 
was a single closed contour, and not the double closed contour seen in Figure 
72. 

To calculate the row of jets in crossflow extensive use was to be made of the 
symmetry boundary condition. The symmetry boundary condition can be applied in 
three basic ways to a row of jets: planes of symmetry placed midway between 
jets on either side of a jet, planes of symmetry placed through the centers of 
an adjacent pair of jets, and, one plane of symmetry placed through a jet with 
another plane placed midway between a pair of jets. For a fixed number of grid 
nodes the latter configuration permits greatest definition of the flow field. 

It was checked that all the symmetry conditions gave the same result. 
Therefore, the most efficient symmetry arrangement can be used without 
problems. This permits maximum jet definition to be achieved for a fixed nodal 
array. Figure 73 shows the symmetry condition used. 
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Fiijure 73 Line ot Jets in a Crossflow - Planes of Symmetry Used in the 
Calculation (lest Case 7 - Khan & Whitelaw) 


Ihe test case flov field was calculated for jets with a pitch to diameter 
ratio of C and a tunnel hei^iht to diameter rat io of 4. Ihe jets had a diaim'ter 
of C.s4 cm and the tunnel had a width of 30 cm. there were 5 holes distributed 
uniformly across the tunnel width, the cross-wind velocity was 10. S ms. The 
calculations were made with a 34 \ 19 \ 11 grid, and were started upstream of 
the line of jets. The hole was represented by 10 \ b grid lines as before. 

f igure 74 presents normal profiles of dimensionless im'>an axial velocities at a 
I. umber of downstream positions both in line with a jet and midway between a 
pair of jets. 
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The results are simultaneously encouraging and discouraging. They are 
encouraging because they are better representations of the measured data than 
the calculations of Khan (Reference 42; and because the profiles for X/D of 4 
to 10, both in line and between jets, agree quantitatively well with 
measurement for Y/D > 1.5. They are discouraging because the inner regions of 
the profile for Y/D < 1.5 for X/D > 4 are not qualitatively correct. 

Detailed consideration of the profiles reveals that close to the tunnel floor, 
both in line with and between jets, the calculated flow is accelerated with 
respect to the cross-wind velocity, and that the flow is not decelerated with 
increasing downstream distance as fast as in the experiment. Deceleration in 
this inner region is calculated to take place, but for X/D greater than 10. 

The mechanism for acceleration of flow near the tunnel floor appears to be the 
flow blockage offered to the mainstream by the jet's penetration and lateral 
spreading which drives the cross-wind toward the wall. This is strongly 
evident at X/D of 4. This diverted cross-wind flow then expands laterally 
behind the jets to maintain the axial velocity there, while itself beginning 
to decelerate, as seen at X/D of 8. 

The axial distribution of grid planes was concentrated about the jet entries 
in order to adequately define the initial events taking place. This results in 
the pleasing profile comparisons for X/D of 4. Downstream, the axial 
distribution of grid planes was sparse. It is therefore believed the delayed 
deceleration of the inner flow arises because of inadequate axial grid planes. 

As this calculation required considerable time to reach a converged solution, 
an economy measure was adopted for subsequent runs. This involved raising the 
criterion for convergence from the 0.5 percent residual source sum normally 
used to a value of 5 percent. 

A comparison of calculated profiles of mean axial velocity at the two levels 
of convergence is enough to justify use of the relaxed criterion for 
subsequent calculations. 

Additional calculations were made for jets with the same pitch to diameter 
ratio of 2, but having a height to diameter ratio of 8. A 34 x 19 x 19 grid 
was used. Profiles of dimensionless mean axial velocity in the plane of a jet 
centerline at a number of axial stations are compared in Figure 75. 

The comparisons of Figure 75 contrast those of Figure 74 in that for the 
height to diameter ratio of 8 the profiles close to the floor are reasonably 
well calculated quantitatively, but the peak of excess velocity is not for X/D 
greater than 4. In the calculations for the height to diameter ratio of 4 the 
profiles close to the roof and the excess velocity peak ware well represented 
while those near the tunnel floor were not. The results are similar to those 
for Cr abb's single jet. 
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Figure 75 (continued) 


Parts of the syminetry boundary condition were checked-out and used in Test 
Case 7. The remaining symmetry condition (upper and lower) was used to 
conserve grid nodes for the opposed jets in cross flow test case, Test Case 8. 
Again, the differences were negligible so the upper and lower symmetry 
condition was also used. 

The case to be considered was for a jet to cross-wind velocity ratio of 2.3 
with the opposed lines of jets 4 jet diameters distant from each other and the 
jets in a line pitched 2 jet diameters apart. Using the symmetry conditions 
the cross sectional calculation area was 2 jet diameters high with a plane of 
symmetry, and its width was 1 jet diameter with planes of symmetry on either 
side. The grid used was 34 x 19 x 11, and the jet origin was represented by 10 
X 5 grid lines, as with Crabb's jet. 

Figure 76 compares measured and calculated transverse profiles of 
dimensionless mean axial velocity at three levels from the floor of the wind 
tunnel, at 6 jet diameters downstream from the hole centerline. 

It is immediately apparent from Figure 76 that although the calculation 
reproduces the shapes of the measured profiles reasonably correctly, there 
appears to be a mass flow discrepancy between the experiment and the 
calculations. 

The cross-sectional area of the calculation domain, which is for one half jet 
and uses symmetry planes, is 2D high by D wide, 2D^. The cross-wind 
volumetric flow into this area is For a velocity ratio uj/U^o of 2.3 

the jet volumetric flow ipto the calculation domain is one half hole area 
7 tDv 8, multiplied by 2.3U*,. Hence, the total volumetric flow is the sum of 
these, or, 2.9032D^U„. The mean velocity leaving the section is therefore 
1.4516Hx,. This simple accounting provides a continuity check on the 
calculations and measurements. The calculated mean velocities leaving the 
section satisfy this limit. Consideration of the measured profiles in Figure 
76 together with the assumption of a smooth function variation between the 
Y-planes, indicate that the experiment does not satisfy this continuity check. 

If the velocity ratio in the experiment were actually unity instead of 2.3, 
then continuity would be satisfied. Kahn (Reference 42} did in fact, make some 
measurements for such a condition, and it is possible some of his figures were 
mislabelled. In view of this problem the case was not proceeded with. 
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Figure 76 (continued) 


8.2.7 Radiation 

Test Case 11 was an exact thermal radiation calculation using the spectral 
method of Semerjian & Segalman. The case used for this study was a reacting 
flow case from Test Case 9 (Roquemore's experiment). The 2D-TEACH radiation 
calculation was compared with the more “exact" spectral calculation made using 
the TIERS code. The radiation fluxes to the wall were compared at two axial 
locations of 25.3 and 74 cm from the face of the flameholder. The tunnel 
operating conditions were 2 kg/s airflow and 16 kg/hr propane. 

Table 15 summarizes the results of the comparison. Two TIERS calculations were 
made: one using the full products for combusting propane in air, and one using 
the "radiation products" assumed in the TEACH radiation calculation, 
consisting of carbon dioxide and water vapor only. 
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Table 15 

Radiation Comparison 


^IQINAL PAGE R5 
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CASE 

F 

L 25.3 cm 

U 

X 74 CIS 
W/Cis^ 



TIERS 

TIERS 

TEACH 

(CO^♦H^O) 

(full combustion products) 

1.489 

0.802 

0.367 

0.090 

0.031 

0.015 


The comparison shows that the neglect of soot radiation together with the 
conclusion that radiation heat transfer only affects the flow field to a minor 
degree for Test Case 9 were justified on a relative basis. However, It does 
suggest that the single flux method used In TEACH Is unacceptably inaccurate. 

The conclusion regarding accuracy Is somewhat surprising In view of the 
previous comparisons that Pratt & Whitney Aircraft had made In Its own 
investigations of the radiation model. These comparisons had shown acceptable 
agreement of the TEACH radiation and gas emissivity models In a flow situation 
very similar to that of Test Case 9, with a Monte Carlo solution. The maximum 
wall heat fluxes agreed within 5 percent, and the emmisivltles with Pratt & 
Whitney Aircraft's empirical gaseous emissivity data base. The overall 
radiation levels in Test Case 9 were very low (less than 1/lOth those of the 
Monte Carlo comparison and less than l/20th those of gas turbine combustors). 
It had been observed in the comparison with the Monte Carlo calculations that 
the single flux model did tend to overestimate at low levels of radiation. 

There remain some questions concerning the setup of the TIERS calculation 
Itself, and also how the TEACH radiation model handles Intervening gas layers. 
More Investigations are required to resolve these questions. 

The TIERS results in themselves indicate that when the overall level of 
gaseous radiation is low (low operating equivalence ratio), the radiation 
calculation should include the full range of combustion products, and not just 
the strong emitters carbon dioxide and water vapor. 

Additional comparisons are clearly needed, particularly at more representative 
radiation levels. 


8.3 Assessment of Model Performance 
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The two equation, K-s turbulence model performs tolerably well for co-axial 
jets in confined sudden expansion and for constant density bluff body flows 
(see References 17, 32 and 80 also). For isothermal variable density systems 
use of a constant turbulent Schmidt number for turbulent mass transport is 
inappropriate (see Reference 93 also). Table 16 gives values for turbulent 
Schmidt numbers that have been found suitable for various flow situations. 
Using hybrid differencing mean concentration profiles can be calculated 
reasonably accurately using recommended turbulent Schmidt numbers provided the 
flow field does not contain a large number of stagnation points. Recirculation 
zone sizes and strengths and mean velocity profiles can be calculated fairly 
accurately. Turbulence quantities are calculated with qualitative accuracy. 


Table 16 

Appropriate Turbulent Schmidt Numbers for Nonreacting Flows 
and Using Hy'»Tiri Differencing and K - e Turbulence Model 


How Type 

Characteristics 

^t 

Boundjry 1«yers 

Predominantly unidirectional 

0.9 

Round free jets 

2-D, external entrainment 

0.7 

Simple recirculation zones 

Simple stagnation systems 

0.5 

Confined round jets 

2-D, self-entrainment 

0.7 

Complex recirculation zones 

2 or 3-D, multiple stagnation 
points, multiple recirculations 

0.2 


While it is difficult to separate all the effects in the swirling flow cases 
studied (see Reference 70), it does seem that the K-e model is inadequate. The 
critical swirl number range for the existence of an internal recirculation is 
not calculated correctly. Excessive viscosity (either from the model or false 
diffusion or both) in the calculation rapidly reduces free vortex swirl to 
forced vortex swirl, see Figure 77, while flows that are naturally forced 
vortex in nature are better calculated. Figure 78. The back-pressure effects 
on the recirculation zone shown in Figure 13 are not calcula^’ed at all. 
Although some swirling flows were calculated tolerably well, the overall model 
performance has to be considered as poor. 
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77 Failure to Calculate Free Vortex/Forced Vortex Core Flows Due to 
Excessive Viscosity 
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Figure 78 Forced Vortex Flows Are Calculated at Least Qualitatively with K - e 
Turbulence Model and Hybrid Differencing 


The liquid fuel spray model has performed well so far as it was tested. 
However, a better statistical approach for the liquid bouiuliry conditions 
appears desirable. 

As anticipated, when compared with a reacting flow of low equivalence ratio at 
ambient operating conditions, the eddy break-up combustion model showed its 
lack of "chemical viscosity" due to the infinitely fast reaction rate 
assumption embodied in it. This is not really so much a fault of the model per 
se, as it is the result of applying it inappropriately. However, it does show 
that for calculations at low engine power levels a model with finite-rate 
chemistry is necessary. 

In the jet- in-crossf low studies using 3D-TEACH, the general features of the 
flows were correct, and some regions of good quantitative agreement were 
achieved. Unfortunately, the calculations were ♦otally dominated with false 
diffusion arising from the coorse grids it was necessary to use, so that 
quantitative accuracy could not really be achieved. However, there was nothing 
in the calculation at this level to i ast doubt on the tuibulence model. 
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The problems of making three-dimensional calculations with the present 
numerics are highlighted in Table 17. On the basis of this table, the code in 
its current form can hardly be considered a viable tool for complete, 
three-dimensional combustion chanter calculations, irrespective of the 
performance or lack of it in the physical modeling. 

Table 17 

Comments on Three-Dimensional Calculations Based on 
Jets-In-Crossf low Experience 


General Observations: 

0 Difficulties in representing circular jet port with rectilinear grid 

0 Convergence rate adversely affected by sudden changes in grid density 

0 Inadequate resolution for even single jet (10,000 nodes, 4 ICG machine) 

0 Excessive CPU times required even with symmetry boundary condition: 

Single jet: 54 mins. (0.5X convergence) 

Line of jets: 111 mins. (PZH4 0 0.5X), 45 mins. (P4H8 0 5 %) 

Opposed lines of jets: 305 mins. (5X convergence) 

0 Excessive numerical diffusion 

Pe > I 2 I everywhere 
Pe > 100 over most of field 
Pe > 1000 in places 


The single-flux radiation model, like the combustion model, was tested using a 
test-case that was somewhat inappropriate and its suitability at more 
representative radiation levels remains to be assessed. It is of course, only 
an interim model. 

Further testing of the combustion, spray and radiation models is indicated as 
necessary. No attempt was made to delineate any effects due to the modeled 
conservation equations being unweighted. This aspect can only be examined when 
some of the other, more important first order, effects have been resolved. 

It should be noted that the plar'^ed step of examining complex flows to 
determine individual model interactions could not be carried out. This was 
because of the problems uncovered with the model performances and the lack of 
suitable experiments. In general, it was not necessary to consider the 
accuracy of the calculations for turbulence quantities as the mean flows 
revealed the problems. Detailed analyses of the reasons behind the physical 
model shortcomings have not been presented. 
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9.0 CONCLUDING REMARKS 


Phase I has proved useful in that it has provided a detailed and systematic 
examination of the physical modeling used, and it has shown the weaknesses 
inherent in this modeling. It has also provided substantial insight into the 
powerful influences exerted on solution accuracy by the inlet boundary 
conditions, and the numerics used in the computer codes. 

The study made it apparent that serious difficulties exist with the present 
generation of computational fluid dynamics codes for application to the 
combustor. These difficulties are not perceived as being fatal however. The 
economic and competitive driving forces in the industry that gave rise to the 
codes are still real and urgent, and the codes have indeed demonstrated great 
potential for productivity improvement in addition to great technical 
capability. Acceptance of computational methods at the working level is 
continuously growing. Ways to improve the existing methods exist, and others 
can be devised. There is no reason or need to abandon the effort. Improvement 
of the codes should be continued. 

It was revealed that the two-dimensional code currently is suitable for 
general engineering use in scale-up, parametric and diagnostic studies. The 
three-dimensional code is best limited to use in skilled hands as a diagnostic 
tool for problems of fairly narrow scope. The lack of a viable design tool is 
a matter of serious concern. Vigorous action is needed to remedy this 
situation. 
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10.0 SUMMARY OF RESULTS 


A comprehensive study was carried out on the physical modeling embodied in 
state-of-the-art computational fluid dynamics codes applied to the gas turbine 
combustor. Although a rather confused and certainly complex picture emerged 
from this study, a number of clear conclusions were drawn, and a summary of 
the results is given below: 

1. A complete set of good benchmark test cases covering the range of 
constituent flows describing the gas turbine combustor is not available. 

2. A hierarchy of physical modeling exists, and the accuracy achieved with a 
given model depends as much on preceding models in this hierarchy as on 
the particular model itself. 

3. The turbulence model is deficient in momentum transport for swirling 
flows, and in mass transport for complicated recirculation zones. 

4. The solution accuracy depends to large measure on boundary condition 
assumptions and to some extent on the calculation grid. 

5. Three-dimensional calculations are totally dominated by the TEACH-numerics 
for practical amounts of computer storage and processing time. 
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NOMENCLATURE 


Of f-oon QuAL; rv 


Commonly Used Symbols in the Text and Figures 


Lower Case Symbols 


flameholder diameter 
mean mass fraction of species 
length scale of turbulence 
mass flow rate of air 


time-mean mass fraction of fuel 


time-mean mass fraction of oxidant 
time-mean mass fraction of products of combustion 
time-mean static pressure 
rac al axis 


axial velocity (along x-axis) 
jet initial velocity 
centerline velocity of fuel jet 
mass average velocity of air 
radial velocity (along r-axis) 
tangential velocity (along 0-axis) 
axis in preferred direction of flow 


droplet drag coefficient 

jet initial diameter 

mean mass fraction of species 

height of tunnel 

heating value of fuel 

specific kinetic energy of turbulence 

latent heat of evaporation for fuel 
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NOfCNCLATURE 


ORIGfNAL PAGE rs 
OF POOR QUALITY 


Commonly Used Symbols in the Text and Figures (continued) 
Upper Case Symbols 

molecular weight 
P pitch of air jets 

time-mean pressure due to turbulence 
R gas constant 

S source term 

SN swirl number 

T time-mean temperature 

Tg gas temperature 

Uoo undisturbed mainstream velocity 

Y vertical distance normal to x-axis 

Z transverse distance normal to x-axis 

Greek Symbols 


r 

8 . 

1 

£ 

0 

V 

P 

‘^t 

n 


turbulent exchange coefficient 
Kronecker delta 
dissipation rate of K 
any scalar quantity 
turbulent or eddy viscosity 
kinematic viscosity 
density 

turbulent Schmidt number 
any independent variable 
void fraction 


overbar denotes time-mean, generally 

' denotes a quantity that fluctuates randomly with time 

Other symbols are defined in the text. 
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EXAMPLE OF RETRIEVAL SYSTEM FOR DATA BASE 


Figure A-1 shows the flow diagrcm of a retrieval system. The following 
comments regarding the classification and contents of the elements comprising 
the data base management retrieval system may be interpreted as applying to 
the description of the files within a virtual machine environment. 

All large computer systems have available powerful data management systems. 
However, these systems often require that the data base be written in a manner 
that is peculiar to a particular machine. The format used for the input 
control and data source elements described below was chosen so as to make the 
data base (viz-a-viz, the computer program to manipulate the base) readily 
usable on other computers. Also, because this retrieval system has been 
designed primarily to be illustrative, little attention has been paid to 
making this version fail-safe, and methods of producing output that can be 
used directly by another computer program have not been included. 



Figure A-1 


Flowcharc Showing Operation of Data Base Management Computer System 
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There are two types of input control elements in the file. The first of these 
element types contains the names of the categories of data contained in the 
data base. The categories currently available are shown in Table A-1. The 
prefix TOC (for Table of Contents) has been appended to each category. Each 
full name (prefix + category) appearing in columns 1 through 12 of the card 
images listed in Table A-1 refers to the name of an element (the second type 
of input control element) to which the program can refer to determine the 
sources of data available in the data base. The retrieval program displays for 
the user the descriptive information contained in columns 21 through 72 during 
the category selection process. An example of the second type of input control 
element is shown in Table A-2. Each entry in this table consists of three card 
images. The first card image contains the name of an element containing all of 
the Information for one source of benchmark quality data. The next two card 
images contain abbreviated bibliographical information in columns 7 through 72 
to aid the user in selecting the appropriate source. The computer program 
displays these latter cards during the selection process. 


Table A-1 

Categories Contained in the Data 3ase 


TOCISOTHRM 

TOCEXOTHRM 

TOCSPRAYS 

TOCKINETICS 

TOCTURBULNT 

TOCSOOT 

TOCRADIATION 


FLUID MECHANICS - ISOTHERMAL FLOW 

FLUID fCCHANICS - EXOTHERMIC FLOW 

SPRAYS, TWO-PHASE FLOWS 

COMBUSTION KINETICS 

TURBULENT COMBUSTION 

SOOT FORMATION AND DESTRUCTION 

RADIATION - ANALYTICAL AND EXPERIMENTAL 


Table A-2 

Example of Table of Contents Element - TOCRADIATION 


RADGODRIDGE 

GODRIDGE.THURLOW, WALLIS - (tTHODS FOR STUDYING FLAME CHARACT. ON 
HEAT TRANSFER IN FURNACES, J.INST.FUEL, 1958 (EXPERIMENTAL)BMARK»3 
RAOGROSSHDLR 

GROSSHANU: ER - RADIATION FROM NONHOMOG. FIRES, FACTORY MUTUAL 
REPT., 1979 (ANALYTICAL),BMARK*8 
RADOSUWAN 

OSUWAN - RADIATIVE HEAT TRANSFER IN CYLINDRICAL TEST FURNACE, 

PHD DISSERTATION, 1971 (EXPERIMENTAL AND ANALHICAL), BMARK«4 
RADSEMERJ 

SEMERJIAN AND SEGALMAN-TURBINE ENGINE INFRARED SIGNATURE (TEIRS) 
PROGRAM, PRATT AND WHITNEY, 1976, (ANALYTICAD,BMARK=9 

RADWU 

WU-COMP. OF PERFORMANCE OF GAS + OIL FLAMES IN CYLINb^ICAL FURNACE, 
J.INST.FUEL, 1969 (EXPERIftNTAL),BMARK=3 


^ f POOR 

If a new source is added to the data base for a particular category, then a 
data source element (described below) is created for the source and added to 
the file and the Table of Contents element for its category is amended. If a 
new category is createo, then the element shown in Table A-1 is amended, the 
data source elements for the new sources are added to the file, and a Table of 
Contents element is created. This simple procedure permits modifications to be 
made to the data without the necessity of altering the computer program. 
Furthermore, portability of the data base is enhanced. 

Data Source Elements 


Each data source element contains the information from a single reference. The 
element consists of card images containing bibliographical and descriptive 
information, commentary and data. Portability of the data source elements is 
assured since card images from one computer system can be readily converted 
for use on another system. 

Each element consists of two major sections. Each card image for the first 
section consists a 6- character keyword followed by up to 66 characters of 
alphanumeric data. The keywords are listed in Table A-3. Columns 73 through 80 
of the standard 80-coiumn card image are not used by the retrieval program 
described below. Keywords less than 6 characters in length are left- justified 
and filled with blanks. Card images containing the same keyword are grouped in 
consecutive order. The keywords can be used to control the flow of the data 
base retrieval program. 


Table A-3 

Keywords Used in the Data Base 


Keyword 

Full-Name 

CATEGO 

CATEGORY 

SUBCAT 

SUBCATEGORY 

SOURCE 

SOURCE 

EXPERI 

EXPERIMENT DESCRIPTORS 

BENCHM 

BENCHMARK EXPERIMENT RATING 

COWCN 

COPWENTARY 

MEASUR 

MEASURED PARAMETERS 

VARIED 

VARIED PARAMETERS 

INITIA 

INITIAL CONDITIONS MEASURED 

INSTRU 

INSTRUMENTATION 

GEOMET 

GEOPETRY 

TOC 

TABLE OF CONTENTS 

ENDKEY 

END KEYWORDS 
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The computer program assumes that the keywords will be encountered In the 
order shown in Table A-3. This restriction is not inherent in the data base. 
Except as noted, there is no restriction concerning the number of card images 
having the same keyword except that these card images must be grouped together 
as noted earlier. Obviously, it is logical to expect that only one card image 
will be encountered for some keywords. A brief description of each of the 
keywords in Table A-3 as interpreted by the computer program will now be given. 

Keywords CATEGO and SUBCAT 


The keyword CATEGO refers to the category, or class, for which the data in the 
source are relevant. For completeness, the data may be further classified by 
subcategory (SUBCAT). The categories shown in Table A-1 in practice refer to 
the subcategories. Only one entry for each of these keywords is expected. 

Keyword SOURCE 

Entries under SOURCE refer to the bibliographical information such as authors, 
source title, archival reference, date of publication, authors' organization, 
etc. 


Keyword EXPERI 


Five entries are made under this keyword. The first entry is interpreted as 
the type of data (experimental, analytical, or both) contained in the 
citation. The second entry refers to the dimensionality of the experiment (or 
analysis). The third entry describes the medium (e.g., air) while the fourth 
card image refers to the injectant (e.g., droplets of JP-4). The fifth entry 
consists of a one-line description of the experiment or analysis. 

Keyword BENCHM 


The single card image using this keyword gives the benchmark rating, 0 through 
10, where 10 represents the highest rating. 

Keyword COMMEN 

This keyword contains commentary that may aid in the interpretation of the 
data or in qualifying the benchmark rating. 

Keywords MEASUR, VARIED, INITIA and INSTRU 


These keywords are used to describe the parameters measured and varied, the 
initial conditions measured, and the types of instrumentation used. 
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Keyword GEQMET 


The geometry used in the experiment or analysis is described in this section. 
In some cases, reference is made to tabular information describing the 
geometry presented under the keyword TOC. 


Keyword TOC 

This keyword refers to the Table of Contents for the source which lists the 
tables and figures available in the data source. The format of each of these 
card images is as follows: 


Cols. 1-6 
Cols. 7-12 
Cols. 13-20 
Cols. 21-72 


TOC ( left- justified) 

TABLE or FIGURE ( left- justified) 

Table or Figure number (integer, right-justif ied) 
Description (alphanumeric data) 


In the present context, a TABLE consists of one or more lines of alphanumeric 
information (e.g., a table of fuel properties); a FIGURE consists of numeric 
data that can be manipulated or plotted. TABLE entries precede FIGURE entries. 
Up to 999 tables and 999 figures may be included in the data source element. 


Keyword ENDKEY 


This keyword is used to terminate the keyword section. 

The second section of the data source element is the data section. The 
information for each table and figure is entered in the same order as the 
tables and figures are listed in the Table of Contents. 

The first line for each TABLE consists of a 72-character alphanumeric title, 
left- justified. The second line contains an integer, M, right-justif ied to the 
tenth column, that indicates how many lines of information follow for this 
table. The next M card images each contain up to 72 characters of alphanumeric 
data, left-justif ied. 

The first line for each FIGURE consists of a 72-character alphanumeric title, 
left-justif ied. The second card image contains the integers, M and N, right- 
justified to columns 10 and 20, respectively. The integer M refers to the 
number of card images that follow the column labels (e.g., the number of 
entries for the independent variable); the integer N refers to the number of 
dependent variables (N< 5). 

The next three card images contain alphanumeric information that may be used 
to label the data columns. Each card image contains N+1 fields of alphanumeric 
data, 10 characters per field. The computer program inserts blanks between the 
fields when printing the data as an aid to readability. Each of the following 
M lines contains a value for the independent variable and the corresponding 
values for the N dependent variables. The format is (6E10.5). 

An example of a data source element is presented in Table A-4. 
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Table A-4 
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Example of Data Source Element 


CATE GO :WDIATI0N 
SUBCAT 

SOURCE GODRIDGE,A.M.,G.G.THURLOW,AND J. WALLIS: 

SOJRCE A METHOD OF STUDYING THE INFLUENCE OF FLAME CHARACTERISTICS 
SOURCE ON HEAT TRANSFER IN FURNACES, J. INST. FUEL. 

SOURCE VOL. 31, NO. 214, NOVEMBER 1958, PP. 491-505, 

SOURCE (BRISTISH COAL UTILIZATION RESEARCH ASSOCIATION) 

EXPERI EXPERIMENTAL 
EXPERI AXISYMNETRIC 
EXPERI AIR 

EXPERI TOWN GAS OR METHANE 

EXPERI RADIATION IN CYLINDRICAL FURNACE WITH CONCENTRIC FUEL/AIR JETS 
BENCHM 3 

COMMEN THE 3-IN.SQ. ANNULAR AIR PASSAGE DUMPS INTO THE 5-IN.DIA. 

COHMEN FURNACE. THE GAS JET DIAMETERS ARE 3/8 AND 3/16 IN. 

COMMEN NEGLECTING THE THICKNESS OF THE GAS INJECTION TUBE, 

COMMEN WHICH PROTRUDES SLIGHTLY INTO THE FURNACE, ONE CAN 
CO^WEN ESTIMATE THE O.D. OF THE AIR PASSAGE. 

COMMEN RESULTS ARE ALSO REPORTED FOR PARTIALLY-PREMIXED MIXTURES 

COMMEN OF FUEL AND AIR. THESE RESULTS ARE NOT INCLUDED IN THE DATA BASE. 
COMMEN A FEW TEST RESULTS WERE OBTAINED COMPARING THE HEAT ABSORBED 
COMMEN PER UNIT WALL AREA FOR TOWN GAS AND METHANE AT NEARLY THE SAME 
COMMEN HEAT INPUT RATE. THE DIFFERENCES DO NOT APPEAR TO BE SIGNIFICANT. 
COMMEN THE WALL HEAT FLUXES ARE DETERMINED FROM THE TEMPERATURE 

COMMEN DROP ACROSS AN ACCURATELY MACHINED SOLID REFRACTORY ROD INSERTED 
COWEN INTO THE WALL. THE EMISSIVITY OF THIS MATERIAL AND THE ROD 
COMMEN SURFACE TEMPERATURE (I.,E.,THE SURFACE OSEENO BY THE COMBUSTION 
COMMEN GASES) ARE NOT REPORTED. THE USEFULNESS OF EITHER THE FLAI€ 

COMMEN RADIATION DATA OR THE WALL HEAT FLUX DATA IN ASSESSTNG THE 
COMICN ACCURACY OF A RADIATION MODEL IS SOMEWHAT LIMITED BY THE LACK 
COMMEN OF INFORMATION ON THE SURFACE TEMPERATURE. 

COItlEN THE AUTHORS ESTIMATE SOOT CONCENTRATION— NOT MEASURED- 

COMMEN FROM THE EMISSIVITY DATA. THEY CONCLUDE THAT THE VALUES OBTAINED 
COHMEN ARE CONSISTENT WITH REPORTED VALUES FOR LOW C/H RATIO FUEL 
COMMEN OILS AND STATE:aTHIS AGREEMENT HELPS JUSTIFY THE RADIATION 
COMMEN MEASUREMENTS IN THIS TRIAL.o FLAME LUMINOSITY IS NEGLECTED 
COMMEN IN ASSESSING THE AXIAL RADIANT HEAT TRANSFER VARIATION. 

COMMEN THE BENCHMARK RATING IS BASED ON (l)THE LACK OF WALL 

COMMEN TEMPERATURE INFORMATION, (2) THE LACK OF INFORMATION 
COMMEN INDICATING WHETHER RADIATION FROM LUMINOUS FLAMES WAS 
COMMEN IMPORTANT. AND (3)THE FACT— COMMON TO MANY EXPERIMENTS— 

COMMEN THAT LOCAL RESULTS MAY BE INFLUENCED BY THE DEGREE OF 
COMMEN MIXEDNESS. 

MEASUR WALL HEAT FLUX 
NEASUR C02 DISTRIBUTION 

MEASUR AXIAL VARIATION OF COMBUSTION TEMPERATURE 
MEASUR PROPERTIES OF IN JECTANT-TOWN GAS (TABLE 1) 

VARIED AIR-FUEL RATIO 
VARIED GAS INJECTOR TUBE DIAMETER 
INITIA FUEL FLO.. RATES 
INITIA AIR-FUFL RATIO 

INSTRU THERMOCOUPLES— 2 PER HEAT FLUX BLOCK 
INSTRU WATER-COOLED HEAT EXCHANGER FOR CALORIMETRY 
INSTRU oMONOa RECORDER AND OFYRITEO APPARATUS FOR C02 
INSTRU SUCTION PYROMETER WITH PT-PT-13tRH THERMOCOUPLE 
INSTRU RADIATION PYROMETERS 

GEOICT CYLINDRICAL FURNACE, LENGTH=40 IN., I.0.=5 IN. 

GEOltT CENTRAL FUEL INJECTOR TUBE, 3/8-IN. AND 3/16- IN. DIA. 

GEOMET CONCENTRIC AIR INJECTOR - 3 SU.IN. (APPROX.) - SEE COMMENTARY 


TOC 

TABLE 

1 

PROPERTIES OF TOWN GAS 

TOC 

TABLE 

2 

OTHER DATA AVAILABLE IN REPORT 

TOC 

FIGURE 

1 

EMISSIVITY VS AXIAL DISTANCE 

TOC 

FIGURE 

2 

WALL HEAT FLUX VS AXIAL DISTANCE. 

TOC 

FIGURE 

3 

FLAIC RADIATION FLUX VS AXIAL DISTANCE 


ENDKEY 
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Table A-4 (continued) 
Example of Data Source Element 


TABLE 1 PROPERTIES OF TOWN GAS 
7 

NET CALORIFIC VALUE (BTU/CU.FT.) DETERMINED FROM CALIBRATION RUNS 

RUN 1 Q=469 

RUN 2 Q=465 

RUN 3 Q=457 

RUN 4 Q=462 

RUN 5 Q=452 

NO OTHER PROPERTY DATA REPORTED 
TABLE 2 OTHER DATA AVAILABLE IN REPORT 
11 
C 

AXIAL DISTRIBUTIONS REPORTED FOR THE FOLLOWING; 
EMISSIVITY-LAMINAR FLOW 

AVERAGE GAS TEMPERATURE— LAMINAR AND TURBULENT FLOW 
RADIATION FROM FLAfE-LAMINAR FLOW 
C02 CONCENTRATION— LAMINAR AND TURBULENT FLOW 
WALL HEAT FLUX-LAMINAR FLOW 
C 

SELECTED INTEGRATED VALUES 

PREMIXED FUEL-AIR RESULTS FOR LAMINAR AND TURBULENT FLOW 
COMPARISONS OF WALL HEAT FLUXES USING TOWN GAS AND METHANE 
FIGURE 1 EMISSIVITY (EPS) VS AXIAL DISTANCE(X) 


A/F 

=5 

=6.5 

=5 

=6.5 

INJ.DIAM. 

=3/16 

=3/16 

= 3/8 

=3/8 

X-INCHES 

EPS 

EPS 

EPS 

EPS 

8.0 

0.0 

0.31 

0.0 

0.0 

13.0 

0.37 

0.33 

0.31 

0.03 

20.5 

0.22 

0.13 

0.44 

0.38 

29.0 

0.44 

0.21 

0.42 

0.30 


FIGURE 2 WALL HEAT FLUX (QWALL,BTU/HR-FT**2) ) VS AXIAL DISTANCE(X) 


A/F 

=5 

=6.5 

=5 

=6.5 

•DJAM. 

=3/16 

=3/16 

=3/8 

=3/8 

;NC> 

QWALL 

QWALL 

QWALL 

QWALL 

0.0 

900.0 

825.0 

500.0 

300.0 

5.0 

1350.0 

1175.0 

950.0 

675.0 

10.0 

2550.0 

2075.0 

1800.0 

1350.0 

12.5 

3100.0 

2500.0 

2250.0 

1725.0 

15.0 

3350.0 

2650.0 

2575.0 

2100.0 

20.0 

3400.0 

3050.0 

2450.0 

2600.0 

25.0 

3275.0 

3275.0 

2300.0 

2675.0 

30.0 

3000.0 

3000.0 

2100.0 

2575.0 

!E 3 FLAME 

: RADIATION FLUX 1 

:qrad,erg/sec-cm**2) 

8 

4 




A/F 

=5 

=6.5 

=5 

=6.5 


INJ.DIAM. =3/16 =3/16 =3/8 =3/8 

X-INCHESqRAD10**-7QRAD10**-7QRAD10**-7QRAD10**-7 
5. 


VS AXIAL DISTANCE(X) 
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Program Operation 

The computer program, control elements, and data source elements are contained 
in a named file. A flowchart summarizing the essential operations of the 
computer program used to retrieve and display the data is shown in Figure A-1. 
The computer program is executed as follows: 

@ASG,A FILE NAME 
0USE FILE, FILE NAME 
0ADD FILE. START 


The element START allocates storage for the logical units (Table A-5) used by 
the program. These logical units include the print file (PRINT) which is 
assigned to logical unit 20; logical units 8 through 14 are scratch 
(temporary) units. The computer program is then started. The element 
describing the categories of data contained in the data base (see Table A-1) 
is then read and stored. 


Table A-5 

Logical Units Used by the Program 


Logical Unit 
Number 

5 

6 
8 

9 

10 
n 
12 

13 

14 
20 


Purpose 

Input unit - terminal 

Output unit ' terminal used for Input 

Data source element contents except keywords 
TOC. COHMEN and ENDKEY and data 

Data source element - COHMEN 

Data source element - Taple of Contents 

Data source element - TABLE data 

Data source element - FIGURE data 

CATEGORIES element (Table A-1) 

TOC category element (e.g.. Table A-Z) 

Print file 


The categories are displayed and the user selects the category that is to be 
examined. The table of contents element for this category is read and stored. 
Then, the program retrieves the information contained in the data base for 
each source listed in the Table of Contents. The first subsection consists of 
the information for the keywords CATEGO, SUBCAT, SOURCE, EXPERI and BENCHM. 
The second subsection lists the entries for the keyword COMMEN. The next 
subsections display the data for keywords MEASUR, VARIED, INITIA, INSTRU and 
G£0^tT, respectively. Finally, the Table of Contents (for TABLE and FIGURE 
data) is displayed. The user then selects the TABLE or FIGURE to be viewed. 
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After displaying each subsection, TABLE or FIGURE, the program pauses so that 
(1) the user may clear the screen if a video terminal without a scrolling 
provision is being used and (2) the user may elect to continue with the 
current source or proceed to the next source. If the user elects to proceed to 
the next source, the program will ask the user if a printed copy of the 
current source is desired. 

After each source is examined, the categories will again be displayed and the 
above sequence may be repeated. When the user terminates the session, the 
program will display the command(s) necessary to obtain a printout and it will 
insert the necessary information into the runstream to detach all logical 
units {except 5 and 6) shown in Table A-5. After program termination, the user 
then logs-orf the computer. 
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